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Abstract 



The so-called Schwinger effect, the creation of electron-positron pairs in strong elec- 
tric fields, has been discussed for nearly eighty years even though an experimental 
verification of this effect seemed to be impossible: The generation of electric fields of 
the order of lO^^V/m, which are required to overcome the exponential suppression 
of the Schwinger effect, was not feasible. Nevertheless there has been substantial 
progress in the theoretical models and methods as well as a rapid development in 
LASER and computer technology during the last decades. As a consequence, plans 
for high-intensity LASER facilities have been made and thus an experimental mea- 
surement of the Schwinger effect could become possible in the near future. 

In this master thesis the non-perturbative electron-positron pair production in 
time-dependent electric fields is investigated. The quantum kinetic formalism is 
employed in order to calculate the electron density for various field configurations. 
The corresponding set of first order, ordinary differential equations is analyzed and 
numerically solved. The focus of this study lies on the dynamically assisted Schwinger 
effect in pulsed electric fields with at least two different time scales. Furthermore, 
interference effects arising in setups with multiple pulses are examined and first 
results for an optimization of the particle number yield by pulse-shaping are given. 
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Introduction 



1.1 Historical Overview 

One of the most important calculations of the 1930s regarding particle production was 
F. Sauters seminal investigation of an electron in a potential of the form V = ux [1] , 
corresponding to a static electric field. His computations within the framework of bar- 
rier scattering revealed that the transmission coefficient is exponentially suppressed 
and becomes sizable only at field strengths of the order of lO^^V/m, 



A few years later W. Heisenberg and H. Euler [2] pubhshed a seminal paper in 
which they investigated the non-linear effect of virtual electron-positron pairs on the 
electromagnetic field as well as electron-positron pair production. One of the results 
worth mentioning is the effective Lagrangian 
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The weak field expansion of this effective Lagrangian yields as leading order 



correction to the Maxwell Lagrangian [3] : 



with a the fine structure constant and Ecru the critical field strength for particle pair 
production 

E„.^^. (1.4) 

The effective Lagrangian (1.3) describes effects, which are not covered by Maxwell 
equations. This lies in the fact that the electromagnetic fields appearing in the 
Lagrangian are of quartic order and thus Maxwell equations in the vacuum become 
modified by non-linear terms. This results in the capability of describing phenomena 
like vacuum polarization or light-light scattering. In particular the amplitude for 
this scattering in the low energy limit is [4] 

In the following years, attempts to formulate a consistent theory of the inter- 
action between electromagnetic fields and matter were made. However, they all 
failed in various aspects. It ultimately took physicists up to the 1950s for a theory 
that could describe the electromagnetic force relativistically and could cope with 
quantum corrections. This theory, whose inventors R. Feynman, J. Schwinger and 
S. Tomanaga received the Physics Nobel Prize in 1965, is referred to as QED or 
quantum electrodynamics, nowadays. 

One of the Nobel laureates, J. Schwinger, also worked on particle pair creation 
in the 1950s. At that time he calculated the probability for particle creation in 
a spatially homogeneous and time-independent electric field using a proper time 
method [5]. This allowed him to derive the probability F of electron-positron pair 
creation per unit time and volume 

F = 2Im(A//), (1.6) 
where the imaginary part of the Lagrangian is given by 

n=l ^ ^ 

The following decades were dominated by the quest for obtaining a standard 
model of particle physics, whereas only minor progress in the field of the Schwinger 
effect was made. Things changed, however, in the 1990s due to the advent of fast 
computers and high-intensity LASERs. Hence within ten years of development strik- 
ing breakthroughs have been achieved in numerical simulations, for instance, for 
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particle production [6-8]. In addition new methods, like the quantum kinetic ap- 
proach [9-12] have been devised based on previous work on quantum field theory 
and transport theory. Recently also experimental physicists have shown an increased 
interest in high-intensity LASER experiments [13] investigating fundamental issues 
hke the Schwinger effect. Most notably, both the XFEL and ELI facilities are sup- 
posed to explore physics at new scales [14, 15]. In this respect, they are going to 
provide the most intense electric fields ever produced in a laboratory in upcoming 
experiments. 

1.2 Aspects of Pair Production 

A direct observation of the Schwinger effect has been considered as impossible for 
a long time as electric field strengths of the order of lO^^V/m seemed out of reach. 
However a modern LASER can provide field strengths of approximately lO^^V/m and 
even stronger LASERs arc planned. Experimental physicists clearly benefit from this 
development as they are provided with a new powerful tool for experiments. The 
challenge from a theorist's perspective is now to perform calculations for realistic field 
parameters. It is obvious, for example, that a LASER field shows a time dependence. 
In addition, the assumption of a spatially homogeneous field simplifies the equations 
a lot, but is also not very realistic. 

The first calculations concerning the pair creation process were formulated as 
a scattering problem. The possibility of an electron tunneling through a potential 
barrier characterizes this effect also nowadays when we talk about the vacuum tun- 
neling probability for an electron. This method for calculation gave a first estimate, 
but soon novel and more reliable approaches have been established. Schwinger for 
example calculated the analytic solution for a time-independent field with the aid 
of a proper time technique, which was for years the most advanced method. More- 
over, semi-classical approaches like the WKB formalism proved to be successful in 
describing the pair production process. In addition for the case of completely uni- 
form electric fields more advanced algebraic approaches have been introduced [16] 
and calculations based on solitons have been made to describe special field configura- 
tions [17,18]. In order to describe more reahstic field configurations, in the 1990s and 
in the beginning of the 21st century new formalisms to calculate the Schwinger effect 
numerically have been introduced, for instance the quantum kinetic approach [11]. 
Other approaches are a closely related scattering-like formalism in terms of the Ric- 
cati equation [19], the Dirac-Heisenberg-Wigner formalism [6,20,21] or the numerical 
worldline formalism [22]. 

Once these methods were established, it became feasible to investigate various 
aspects of pair creation for the first time. A large and growing body of literature has 
investigated time-dependent electric fields. One prominent example is the sinusoidal 
field E = eEcritsiia it/r) [23-25], which was first discussed in the 1970s [26]. Another 
field configuration widely explored in various computations is the Sauter field 
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E = eEcritSech^ {t/r) [27,28], which was also first investigated in the 1970s [29] 
and further analyzed in an S-matrix formalism [30] as well as in the worldline ap- 
proach [22,31]. Moreover, the question has been raised whether it is possible to lower 
the threshold for particle production in order to abet the experimental verification 
of the Schwinger effect. Therefore it has been suggested [28,32-35] that the parti- 
cle production rate could be significantly enhanced by the superposition of electric 
fields with different time scales. Moreover there have been first investigations of the 
Schwinger effect in space-time dependent electric fields [36-38]. Furthermore there 
was theoretical progress as well: For instance, it was shown [39] that there is a dif- 
ference between particle production rate and the rate associated with the decay of 
the vacuum. 

As already mentioned, the WKB approach has turned out to be a very successful 
tool for describing pair production since its first application to the Schwinger effect. It 
is still widely used and gives new insight into this process [19,40,41]. For example, the 
interference effects occurring in pulsed fields with subcycle structure, first observed in 
F. Hebenstreit et al. [42], could be explained conclusively within the WKB approach 
[43,44]. It is also important to state here, that the different approaches give similar 
results [45] and are sometimes even equivalent [19,36,46]. 

The interested reader is referred to the articles [4, 14, 15,47,48] for further infor- 
mation on the topic. If one searches for an extensive review on particle production, 
high-intensity LASERs and its apphcations we encourage the reader to have a look 
at [3,49]. 



1.3 Objective of this Thesis 

We will entirely focus on the quantum kinetic approach in this thesis, as it provides 
a flexible and well-established formahsm. This thesis is divided into three parts. In 
the first part, we introduce our notation and state our assumptions. Afterwards we 
derive a quantum kinetic equation for the particle distribution function 

F {q, t) = W {q, t) f dt'W {q, t') {1 - F {q, t')) cos {26 {q, t, t')) , (1.8) 

with W {q, t) , 9 {q, t, t') as input functions. Due to the fact, that the formulation in 
terms of this integro-differential equation is cumbersome in numerical simulations, we 
will also derive the equivalent set of ordinary differential equations of first order [9] 
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where G, H are auxiliary quantities and u input functions. Then we will calculate 
an analytical solution for the electric field 

E ^ eEcritsech^ {t/r) . (1.10) 

The chapter ends with a presentation of an optimization algorithm. 

Throughout the second and third part we shall work with the field configuration 

A{t) = J2 -Sinte^nh (^^^^ , E{t) = ^e.sech^ (^^) • (l-^) 

The second part assesses different numerical methods, such that one can use the 
results presented as guideline data. 

In the final part all results are shown in a structured way. Firstly, we try to 
find new enhancement effects resulting in threshold lowering. To do that we give 
an account of our attempts to combine different time scales to obtain optimal pulse 
shapes. Second, we point out that one can observe interferences in the particle 
distribution depending on the diverse input parameters. Consequently we are able 
to perform calculations, which show the possibility for an indirect measurement of 
the Schwinger effect. To conclude, the results obtained from applying optimization 
procedures to our system are presented. 
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Quantum Kinetic Equation 



Various methods have been developed to calculate the particle distribution function 
for the created electrons until today. In this chapter the quantum kinetic formalism 
is presented, because it has a number of attractive features. First and foremost it is 
highly flexible and easy to implement in computer programs. Secondly, it is a fast 
way to obtain results for different setups in the parameter region we are interested 
in. 

The derivation of the quantum kinetic equation follows closely the PhD-thesis of 
F. Hebenstreit [37]. It can also be found in the publications [11,27,28]. However, 
we do not employ their compact form, but write it down in a more elaborate way. 
Hence we will explicitly examine the derivation of the quantum Vlasov equation and 
subsequently an analytical solution for the Sautcr field configuration is presented. At 
the end of this chapter we shall also give an introduction to optimal control theory 
and its applications to pair production. 



We start with introductory remarks on the notation for the purpose of a clear struc- 
ture and better readability of the mathematical intense part of this thesis. 

We use the common notation for indices. This means that for Latin letters 
i,j = 1,2,3 hold, while for Greek letters /j,,!/ = 0,1,2,3 is assumed. Moreover we 
use Einstein summation convention 



2.1 Notation 



3 




(2.1) 



/i=0 



to improve readability. 
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In the following we will define basis vectors and introduce the gamma matrices 
7^. We specify the unit vectors 



ei 



/1\ 




\V 



62 



1 





/0\ 


1 

\V 



64 






VJ 



(2.2) 



and denote the unit matrix as I. At certain points we write 1^ to denote the unit 
matrix in a specific dimension. Moreover, we choose the Weyl or chiral basis [50] for 
gamma matrices: 



I2 
I2 



r 



-a' 
a' 



I2 
-I2 



(2.3) 



with the Pauli matrices cr* 

1 

1 



-i 

1 



1 
-1 



(2.4) 



The first foTir gamma matrices form a Clifford algebra, satisfying the important 
anti-commutator relation 

with the metric tensor 



^'^'^ = diag (1,-1, -1,-1). 
Accordingly, they are hermitian and anti-hermitian, respectively 

[1 ) =1 , [1)^-7, (7 j = 7 

and fulfill 



(2.5) 
(2.6) 



(7°) =1, {Y)=-l (7^)' 



I. 



We denote 



rr'f^ ( ™3\ 

and introduce the shorthand notation 

d „ d 



(2.7) 

(2.8) 
(2.9) 

(2.10) 



To distinguish scalars from vectors we write an extra ^ above a vector. For 
adjoint variables we use the notation ^, which is equal to = ^t^o _ (^*)^^o^ 
where means transpose and * complex conjugate. 

The preferred system of units in high-energy physics is the Heaviside Lorentz 
system. Additionally we use natural units and set h — c — 1. In chapter 3 we also 
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set rrie — 1 and express all variables in terms of the electron mass. The electron 

charge e in this system of units is given by p = \/A7ra ~ 0.3028, where a is the fine 
structure constant a ~ 1/137. Last but not least the Maxwell equations take the 
form 

V-E^p, (2.11a) 
V-S = 0, (2.11b) 

V X ^ + (dtB^ = 0, (2.11c) 

V xB - (dtE'^ ^j. (2.11d) 

It is not possible to solve the problem of electron-positron pair creation in the 
framework of QED for any given field configuration exactly. As a consequence, we 
have to employ several approximations: In recent studies theorists concentrated on 
spatially homogeneous time-dependent fields and so do we. If we intended to abandon 
the assumption of spatial homogeneity, we would have to use a more general approach 
such as the DHW-formalism [6,37]. As it is our aim to describe processes occurring at 
field strengths of about Ecr — l.S-lO^^V/m [2] a mean field description of the electric 
field should be a good approximation. In further agreement with the literature we 
treat the electric field classically and quantize only the matter field. 

Motivated by standing waves, we will consider only electric fields whereas mag- 
netic fields are totally neglected 

— * 

S = Vxl=0, ^=-V0-^, (2.12) 

ot 

— * 

with A the vector potential and the scalar potential. Then we exploit the gauge 
freedom 

A^^A^ + d^p (2.13) 
to choose Weyl or temporal gauge = giving us 

l^l+Vp = A 0^0-^ = 0. (2.14) 

We assume the electric field points into the same direction all the time. Thus 
potential and electric field take on the form 

A{t) = (0,l(t)) = (0,0,0, ^,(t)), (2.15a) 
E (t) = (O, E (i)) = (0, 0, 0, E (t)) . (2.15b) 
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2.2 Quantum Vlasov Equation 



As we will see in the following, the description of the whole pair creation process is 
similar to characterizing a many-body problem by a Boltzmann equation 

dtF{q,t)^S{q,t)+C{q,t). (2.16) 

On the left-hand side the temporal change of the distribution of the particles is 
written. On the right-hand side we have both a source term S {q, t) for particle 
production and another one that describes the collisions C {q, t). 

The derivation of the source term S {q, t) is our main goal in the next sections, 
while we will disregard the collision term. This neglect can be justified by applying 
a simple analysis based on the relaxation time method [9] . The idea is to model the 
collision term as 

Cm=^^^MA^^, (2.17) 

TV 

with Feq {q, t) the equilibrium distribution function for fermions. The relaxation time 
Tr is interpreted as the average time between two particle collisions. Due to the fact 
that the density of the created particles is sufficiently low for subcritical electric 
fields, one finds ^ 1 and, as a consequence, C{q,t) <C S {q,t). The dominance 
of the source term over the collision term legitimizes the neglect of collisions in the 
parameter region under consideration. 

2.3 Equation of Motion 

In this section we derive the quantum Vlasov equation. As we do not give an intro- 
duction to quantum field theory, we refer the interested reader for details to books 
on this subject [50-53]. 

A good starting point for the derivation is the QED Lagrangian 

C^i^ (7^9^ _ ie^M^) - - ^F^uF''^ (2.18) 

with charged spinor fields ^(,x,t), ^{x,t) = ^t^o^ rpj^g 

mass term m explicitly 

denotes rriel. Our first objective is to find an equation of motion for the fields ^ (x, t) 
and ^ {x,t). This is achieved by applying a variational principle to the Lagrangian. 
This gives the famous Euler-Lagrange equation of motion 

For the adjoint spinor field this leads to 

(1^-7^) - e^-ff^Ai, + m^' = 0. (2.20) 
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Complex conjugation yields for the spinor field 

17^^(9^^ + - = 0, (2.21) 

(i-f'^dt + i7 • V + e-f^A^ - m) * = 0. (2.22) 

Now we use the specific form of the vector potential in (2.15) and obtain an equation 
of motion for ^ (x, t) 

(^i7°9t + if • (^V - ieA(t)^ - rnj ^ (x, t) = 0. (2.23) 

The spinor fields ^ {x, t) and ^ (x, t) are treated as dynamically independent fields. 
For the sake of completeness we can write down the canonical momenta 

n = ^ = i*t (2.24) 

and the equal time anti-commutation relation (ETAR) 

{^(f,t),n(f,t)} = i5(f-f). (2.25) 

Because of the spatial homogeneity of the vector potential, it is indicated to decom- 
pose both the spinor field ^ (x, t) and its conjugate field 11 (if, t) into their Fourier 
modes 

* = / ^) ' n (f , t) = y ^e-'^-TT (g, t) . (2.26) 

Similar to (2.24) we have for the mode functions 7r{q,t) — itp^ {q,t). The mode 
functions also fulfill an ETAR 

{V^ (q, t) , TT i^, t)} = i {27rf S (q- q') . (2.27) 

Inserting the Fourier expansion into (2.23) and using the relations 

i7° {dt^ (f , t)) = i7° / T^e^"^-^ (g, t)) , (2.28a) 

i (7 • V,) ^ (f , t) = i7 y (ig) e^^-> (g, t) , (2.28b) 



we obtain 



^e'^-^ (i7° (9,V (g, *)) - 7 • (? - el(i)) V - {q, t)) = 0. (2.29) 



(27r) 
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This equation is not only true for the whole integration but for every integrand. 
Hence it is sufficient to work with 

(i7°ai - 7 • (q- el(t)) - m) ^ {q, t) = 0. (2.30) 

Introducing the kinetic momentum '7T{q,t) — q — eA (t) and especially 773 {q, t) — 
qs — eAz (t) the equation of motion for the Fourier modes reads 

(i7°at - 7 ■ 7f {q, t)-m)i{j (g, t) = 0. (2.31) 

Apparently this is an ordinary differential equation of first order. In order to solve 
it, we introduce a quasi-particle representation in the following. 

2.3.1 Quasi- Particle Representation 

The first step in transforming equation (2.31) is to introduce an ansatz for the mode 
function 

V (g, t) = (i7°at - 7 ■ vT (q, t) + m) $ (g, t) . (2.32) 
The equation of motion then takes the form 

(i7°a, - 7 • ^ (q, t) - m) (i7°a* - 7 • ^ (?, t)+m)^ {q, t) = 0. (2.33) 

Exphcit calculation of the left hand side gives 

{—Wf — 17 • ■K^'^dt — im^^dt — i'f'^dt (7 • tt) 

+ (7 • 7f)^ + • 7f + im7°at - • tt - m^) $ (q, t) = 0. (2.34) 

We see that some of the terms cancel each other upon using the relations 

dtil-^ {q, t)) = 1% {qs - eA, (t)) = e^^E (t) , (2.35) 

(7 • Tf {q, t)f = -Inl (q) - Itt^ (q) - IttI {q, t) = -In' {q, t) . (2.36) 
Together with^ 

- (7 • n (q, t)) $ (q, t) = i (7 • ^ (q, t)) 7° (9,$ {q, t)) + iej^E (t) $ (q, t) (2.37) 
equation (2.34) is simplified to 

{W^ + ie7°7^£; (t) + In' {q, t) + m') $ {q, t) = 0. (2.38) 

1 -ij°dt (7 • 7f {q, t)) $ {q, t) = -i7° (7 • 7f {q, t)) (5*$ (g, t)) - i7° {3^ (7 • 7? % i))) $ {q, t) 
= i (7 • % t)) 7° (g, i)) + ie7°7'^ (t) $ (9, *) 
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Introducing the notation uP" {q, t) — tt^ {q, t) + ml this equation is written as 

{W^ + ie7%3E (t) + luj^ {q, t)) $ {q, t) = 0, (2.39) 

which is formal equivalent to an oscillator-type equation with an additional time- 
dependent term. Now wc expand the function $ [q, t) in a basis of eigenvectors of 
^0^3 rjj^g matrix representation of ^^^7^ is 



7%^ 



I \ / 

1 



-a 





-(73 ' ' 



(2.40) 



which makes it easy to determine the eigenvectors. They are given by 



Ri 



n\ 






i?2 = 






VV 



1 




Ra — 





1 



(2.41) 



where Ri, R2 have eigenvalue A = 1 and R3, R4 have eigenvalue A = — 1. This gives 
us the possibility to expand $ {q, t) in the basis of the eigenvectors 



(2.42) 



r=l 



with Xr{Q,t) as time-dependent coefficients^. Inserting this expansion into (2.39) 
leads to 



{Id^ + ILU^ {q, t) + ie7%3E (t)) J2 Xr (q, t) Rr 



r=l 



= {W^ + Iw' {q, t))J2xr {q, t) Rr + ieE (t) ^ Xr (q, t) {^^Rrl . (2.43) 



r=l 



r=l 



The term (*) is an eigenvalue equation yielding Ri^2 and — -^3,4, respectively. We ob- 
tain one equation for every component, which leads to two different scalar equations 



{d^ + {q, t) + ieE (t)) Xr {q, t) ^ r = 1, 2, (2.44a) 
{d^ + (g, t) - ieE (t)) Xr (q, t) ^ r = 3, 4. (2.44b) 

Each of these equations appears twice. Note that the Dirac equation (2.31) is a 
first order differential equation having four linear independent solutions. On the 

^ We will refer to Xr {q, t) as mode functions from this point on. 
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other hand, (2.39) is a second order differential equation such that we obtain eight 
hnear independent solutions in total. However, this redundancy can be removed by 
choosing only one complete set of solutions corresponding to either (2.44a) or (2.44b). 

Hence we decided to skip the second equation. This leaves us with one differ- 
ential equation of second order with fundamental system {q, t) and coefficients 

(q) , c~ (q) where r = 1, 2. Now we have four different components and can write 

2 

^ = E (''r (q) (q, t) Rr + c; (q) x; {q, t) Rr) . (2.45) 

r=l 

Note that the summation goes from one to two and not up to four anymore. Using 
the equation above one obtains 

{q, t) = -7-7? {q, t) + m) $ {q, t) 

2 

= (i7°5i - 7 ■ 7f {q, t) + m) (c+ {q) x+ {q, t) Rr + c; (g) x; {q, t) Rr) (2.46) 

r=l 

and thus the spinor field takes the form 

2 

^ i<l^ t) = Yl K ^) io) + Vr i-q, t) c; {q)) (2.47) 

r=l 

with 

Ur {q, t) = (1709^ -^■^{q,t)+m) [xt t) Rr) , (2.48a) 

Vr i-q, t) = (17°^^ -j-7f{q,t)+m) {x; {q, t) Rr) . (2.48b) 

Moreover, we redefine the coefficients 

ar{q)^4{q), bl{-q) ^ c; {q) , (2.49) 

such that we can write 

2 

^ (?, t)^Yl K t) (q) + Vr i-q, t) bl i-q)) . (2.50) 

In the equation above Ur{q,t) and Vr{—q,t) are time-dependent spinors, whereas 
Or (q) and 6j (—q) can be interpreted as annihilation and creation operators. These 
creation/annihilation operators are fermionic operators and thus the following ETAR 
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hold 

{ar {q) , al (g^) } = {w {q) , bl (g^) } = {2nf [q- q') , (2.51a) 

[ar (q) , as {q) } = (q) , h [q') } = 0, (2.51b) 

{4 (q) , at (g^) } = [bl {q) , bl {^q') } = 0. (2.51c) 

Investigating the Hamiltonian of the system [28] we observe non-vanishing off di- 
agonal elements, which account for particle creation and annihilation. In order to 
calculate the spectrum, we have to diagonalize the Hamiltonian. This is achieved by 
a basis transformation, where we introduce new time-dependent operators Ar {q, t) 
and Br {—q,t), which also show fermionic behaviour 

[a {q, t) , Al i^q', t)]^[Br {q, t) , Bl {^q', t) } = {2nf SrsS {q- q') , (2.52a) 
[Ar (q, t) , As (q', t) } = [Sr (q, t) , (g^, t)] = 0, (2.52b) 
{4 (g, t) , Al (g^, t)] = {Bl {q, t) , Bl (g"^, t) } = 0. (2.52c) 

The relation between the operators (g) , bl (g) and Ar (g, t) , B^ (— g, t) is given by 
a Bogoliubov transformation 

Ar (g-, t) = a {q, t) ar {q) - 13* {q, t) bl {-q) , (2.53a) 
Bl {-q, t)^/3 {q, t) ar (q) + a* {q, t) bl {-q) . (2.53b) 

Note that the coefficients a (g, t) and /3 (g, t) in the equation above satisfy the relation 

\a{q,t)\'+\P{q,t)\' = l. (2.54) 

We will refer to the new basis as "quasi-particle representation" and to the old one 
as "particle representation" from here on. This notion is introduced, because for 
vanishing electric fields the operators a (g) and b^ (g) correspond to single particles, 
in the sense of asymptotic states. On the other hand, the operators A (g, t) and 
B{—q,t) are mixtures of the original operators a{q) and b^ (q)- However, we call 
them quasi-particle operators due to the fact that the quanta corresponding to these 
operators behave like particles with dispersion relation u {q,t). 

As the Bogoliubov transformation gives us the desired change of basis, we may 
express ip (g, t) in terms of new operators as well 

2 

iP {q, t) = Y, {Ur {q, t) Ar {q, t) + K (-g-, t) Bl {-q, t)) . (2.55) 



r=l 



The spinors in the quasi-particle basis Ur (g, t) , K {~q, t) are denoted by upper case 
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letters to distinguish them from the spinors Uj. {q, t) and Vr {—q, t). They are given 

by 

Ur (g, t) = (7°a; (g, t) - 7 ■ 7f (g, t) + m) («:+ (g, t) Rr) , (2.56a) 
Vr i-q, t) = (-7°cc; (g-, t) - 7 • tt (g-, t) + m) {k; {q, t) Rr) . (2.56b) 

The functions nf arc the mode functions in the quasi-particle representation, which 
are chosen according to the ansatz 

(g, t) = =, (2.57a) 

^y2uJ (g, t) {UJ (g, t) - TTg (g, t)) 

K- (g, t) = =, (2.57b) 

^/2u{q,t) {uj{q,t)+m{q,t)y 

with the dynamical phase 9 (g, t) = J^* a; (g, r) dr. Note that cu (g, r) is time-dependent 
due to its dependence on the vector potential A (t) . 

The two functions chosen such that they coincide with the mode functions 

Xo in the case of a vanishing vector potential. The calculation of Xq , which is done 
in A.l yields 

Xo (g", t) = , =e^^"(^"^*. (2.58) 

V2cj {q) {uj {q) ^^3) 

In the case of vanishing electric field and vector potential the following simplification 
is obtained 

u{q,t)^u{q)^ml + f. (2.59) 
Hence the integration in the dynamical phase {q,t) yields 

e{q,t)^ [ uj{q,r)dr^uj{^{t-to). (2.60) 

J to 

We skip an arbitrary phase e'^^^*° and obtain 



e 



4 = xo^ (g', t) = , ,^ =, (2.61) 

^y2uJ (g) {uj (g) ^ gg) 

proofing that our ansatz for {q, t) is valid. 

To summarize, we can express the spinor -0 (g, t) in both the particle and in the 
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quasi-particle representation 



{q, t)^J2 K (^"'' + (-^"'' ^) (-^) (2.62) 

r=l 

2 

= i^r {q, t) Ar {q, t) + Vr i-q, t) Bl {-q, t)) . (2.63) 



r=l 



In order to fully specify the Bogoliubov transformation, we still have to derive an 
equation for a{q,t) and f3{q,t). However a calculation of these two coefficients is 
rather challenging but done explicitly in A. 2. Here it is sufficient to know that these 
coefficients are given by 

a {q, t) = i^u^ {q,t)-7Tl{q,t)ti~ {q, t) {dt - icj {q, t)) x+ {q, t) , (2.64a) 
P (q, t) = -i^ujHq,t)-7ri(q,t)K-^ (q, t) (dt + iu (q, t)) x+ {q, t) . (2.64b) 

An important aspect of the mode functions in the two different representations 
is their similar behaviour for vanishing potentials and fields at asymptotic times 
t — )■ — oo. In this limit, both types of mode functions {q, t) and {q, t) coincide 
with the free mode function Xo {Qj ^) to an arbitrary phase. To check this we may 
calculate the Bogoliubov coefficients in this limit and find 



a{q,t^ -oo) = i\Juj^- g|Xo (^t " M Xo 



iX^ (-2ia;) -— i== = 1 (2.65) 
" 2u^u^ - qi 



and 

|2 1 I , m2 



1/3 {q,t^ -oo)r ^l-\a{q,t^ -oo)r = 0. (2.66) 
This shows that the two basis indeed coincide at asymptotic times i — > — oo. 



2.3.2 Deriving the Source Term 

This is now the final part of the derivation, namely the calculation of the equation 
of motion for the one-particle distribution function F {q,t). This function gives the 
quasi-particle number density for any canonical momentum q. It is a time-dependent 
quantity and defined as the expectation value of the quasi-particle number operator^: 

F (g, t) = ^lim - ^ (4 (g, t) A, {q, t)) . (2.67) 

^ In other words the one-particle distribution function holds as an expectation value of the quasi- 
particles in our system divided by its volume. This is due to the spatial homogeneity of the system. 
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We start with the definition of the particle number density in the time-independent 
basis 

Nr (q) = ^lim i (aj (q) a, (q)) = ^Im i (bl {-q) K {-q)) . (2.68) 

Using the operators Ar {q, t) and Al (g, t), which have been defined in (2.53a) 

Ar (q, t)=a (q, t) a, (q) - /3* {q, t) h]. {-q) , (2.69a) 
4 {q, t) = a* {q, t) al (q) - /3 {q, t) K {-q) . (2.69b) 

the computation of their product gives 

AlAr = \af alar - ^abrar - a*P*albl + \pf brb^ (2.70) 

Correspondingly the expectation value is given by 

(4A> = \af (alar) + \/3f {b^bl) . (2.71) 

As the order of the operators is important we have to rearrange the terms in the 
equation above. For that reason we use the anti-commutator relation defined in 
(2.51a) to get 

{brbl) = {brbl + blbr - ^Jfe.) = bl}) - {b^,br) ^V- {b^) . (2.72) 

We proceed with the calculation 

i {AlAr) = M! l^ala.) + M! (6,6^) = \af AT, + M! (v^ _ ^tX)) 

= (1 - \f3f) Nr + \pf-^-^ (blbr) = Nr+ (1 - 2Nr) , (2.73) 
leading to the distribution function 

2 

F {q, t)^Yl (q) + \P {q, t)f (1 - 2Nr (q))) . (2.74) 

The time-independent quantity Nr (q) specifies the initial conditions of our sys- 
tem. We choose Nr{q) = to get pure vacuum for t — >■ — oo. Thus the one-particle 
distribution function takes the form 

2 

F {q, t) = J2\P {q, t)f = 2\(3 {q, t)f . (2.75) 

r=l 
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Another useful quantity is the quasi-particle correlation function 



^ *^ = i^oovt^ (^1 (-^^ 4 ^)> ■ (2.76) 



r=l 



Again we reformulate this equation to express C {q,t) in terms of the coefficients 
a {q, t) and /? {q, t). This time we work with the operators 

Al {q, t) = a* {q, t) a\ (q) - (5 {q, t) 6, (-g) , (2.77a) 
Bl i-q, t)=f3 (q, t) ar (q) + a* (q, t) bl (-q) (2.77b) 

and multiply them to obtain 

BlAl = a*/3aral + a*a*blal - /5/3a^6^ - a*f3blbr. (2.78) 
Then we take the expectation value 

(BlAl) = Pa* (aral) - a* (3 (bX) . (2.79) 
Reordering of the operators and division by the volume gives 

^ (BlAl) = Pa* (^1 - i (alar)^ - ^ {bX) . (2.80) 

Taking the limit on both sides of the equation results in 

lim (BlAl) = Pa* - - a*pNr. (2.81) 

V^oo V 

Finally we use the initial condition Nr{q) — again, which leads to 
C (g, t) = hm 1 Yl (^l (-^^ ^) 



r=l 



= J] /3 {q, t) a* {q, t) = 2a* {q, t) /3 {q, t) . (2.82) 



r=l 



A formulation of F {q, t) , C {q, t) in terms of a {q, t) and /3 {q, t) is not very ad- 
vantageous as one would have to know the exact mode functions x"*" {q, t). Therefore, 
we try to obtain an equation for F {q, t), which is easier to compute. After another 
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tedious calculation, performed in A. 3, we obtain 

« iQ, t) = {q, t) 15 {q, t) e^^^®*), (2.83a) 

$ (q, t) = -lw (q, t) a {q, t) e-^^^(«"*'*), (2.83b) 

with 

W {q, t) = ^^f^t;^^^^ , {q, t)=ml + {q, t) , (2.84) 

eliq)^ml + ql + ql (2.85) 
Taking into account these equations, we find^ 

C (q, t) = -W (q, t)(l-F (q, t)) e-^^^(^'*) (2.86) 

and^ 

F {q, t) = -W {q, t) Rc {C {q, t) e^''^'^''^) . (2.87) 
Integration of C {q,t) with respect to t gives an expression we can insert into (2.87) 

F {q, t) = -W {q, t) Re (C {q, t) e^^^^^'*)) 

= W {q, t)Re(^J dt'W {q, t') {I - F {q, t')) e^iC^C^*)-^®*'))^ . (2.88) 

Introducing the shorthand notation 

e{q,t) -e{q,t') = [ oo{q,t")dt" = e{q,t,t') . (2.89) 
Jt' 

we finally arrive at 

F (q, t)^W (q, t) f dt'W (q, t') (1 - F (q, t')) cos {26 (q, t, t')) . (2.90) 

Given the one-particle distribution function, the asymptotic particle number per unit 
volume can be calculated according to 



/oo 
d^qF{q, 00). (2.91) 
-00 



1 (7 = 2 (a*p + a*/?) = 2 (11^/3*6-2^^ - la*Wae-'^''^) = -W {a*a - /3*/3) e-^^^ 

= -W - 1^1^) e-2'« = -W^ (1 - 2 \l3f^ e'^'^ = -W {1 - F) e'^'^ 

2 F = 2dt = 4Re (^^*) = 4Re {-ll3Wa*e^'''^) = -WRe (2a* ^e^'^) = -WRe (Ce^^^) 
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Equation (2.90) is an integro-differential equation, where the right hand side 
acts as a source term for the system. The factor (1 — F {q,t')) under the integral 
gives evidence for non-Markovian behaviour. Furthermore, this term represents Pauh 
exclusion principle, thus we refer to this term as Pauli-blocking factor. The second 
interesting part is the cos (26' {q,t,t')), which is an indicator for a highly non-local 
problem in time. 

We have not taken into account the electric fields, which are created by the 
produced particles, so far. This disregard is based upon previous investigations [9,27], 
where it has been shown, that those internal fields only play a role for critical field 
strengths or high particle densities. However, their contribution to the final particle 
density is rather low even in these regimes. Therefore we decided to entirely neglect 
this back reaction effect in this thesis. 



2.3.3 Differential Equation 

In order to solve the integro-differential equation (2.90) we have to calculate the 
occurring integrals as exact as possible. In this respect the highly oscillating term 
cos {26 {q,t,t')) is numerically challenging. Moreover, due to the non-Markovian 
behaviour, the complete time-history of the one-particle distribution function has to 
be taken into account, which is potentially error-prone. To avoid these complications 
we rewrite this equation as a coupled differential equation as proposed by J. C. R. 
Bloch et al. [9]. In close analogy to this paper we introduce the auxiliary functions 



G{q,t)= f dt'W {q,t') {1- F{q,t')) cos {29 {q,t,t')), (2.92a) 

''Wac 

H{q,t)= f dt'W{q,t'){l-F{q,t'))sm{2e{q,t,t')), (2.92b) 



where 



4(^)=m^+??+?2^ (2.93) 

e {q, t, t') = ^* CO {q, t") dt", {q, t) ^ ml + (q- el(i))' . (2.94) 

The derivatives of the auxiliary functions with respect to t yield 

G {q, t)^W {q, t){l-F {q, t)) - 2u {q, t) H {q, t) , (2.95a) 

H{q,t)^2u{q,t)G{q,t), (2.95b) 

where we have used 

e {q, t, t') = dt (^J^^ CO {q, t") dt"^ = CO {q, t) . (2.96) 
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Combining the equations (2.95b) (2.95a) with (2.90) gives us now the possibihty to 
work with a coupled differential equation of first order 

F{q,t) = W{q,t)-G{q,t), (2.97) 
G (g, t) = W (q, t) (1 - F (g, t)) - 2uj (g, t) ■ H (g, t) , (2.98) 
H{q,t) = 2uj{q,t)G{q,t). (2.99) 

Its matrix form is 




(2.100) 



which indicates a special structure of the system. Initial conditions are chosen such 
that one starts with a pure vacuum at asymptotic time t — >■ — oo 

F{q,t^ -oo) ^G{q,t^ -oo) ^ H {q,t ^ -oo) = 0. (2.101) 



2.4 Analytical Solutions 

The relationship between time-dependent electric fields and the particle number den- 
sity has been widely investigated in the last decade. One of the most-studied field 
configurations is the pulsed Sauter field 



A (t) = -eEcrTtanh j , E (t) = eEcrSech^ j , (2.102) 

as it allows for an analytical solution. As nearly all of our investigations in chapter 
3 and chapter 4 will be based on combinations of such simple pulses, we give the 
complete calculation of the one-particle distribution function F {q, t) for a single pulse 
in this section. 

First we set rrie = 1 and introduce the notation 

A {t) = eA (t) = -ertanh (^^^ , (2.103) 

E (t) = eE (t) = £sech^ (^^^ " ^ ~ (r ) ) ' (2-104) 

For convenience we redefine A{t) — A (t) and E (t) — E (t) to keep the notation of 
the variables simple. This redefinition holds for the rest of this section, giving 

{q, t) = (^3 - A, {t)f , cj^ {q, t) = 1 + qf + qj + tt^ {q, t) . (2.105) 
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Consequently we obtain for the equation of motion (2.44a) 

(..= ..^ff,)+.sech^(i)).(,,)^0. (2.106) 

where we have skipped the irrelevant index of the mode function x {q, t)- For conve- 
nience, we transform the time variable 

t ^ M = J (1 + tanh ( - ) ) , [-00, oo] ^ [0, 1] (2.107) 



2 V V'^. . 

which leads to 

A{u)^-eT{2u-l), E {u) ^ 4:ue {1 - u) . (2.108) 
Additionally, we calculate 

t = rarctanh {2u - 1) , (2.109) 
lA •2ft\\ 1 u2^t 



u{l-u) = -{1- tanh^ ( - ) ) = ysech^ ( - ) , (2.110) 



dt = {dtu) du = ^sech^ (-] ^-u{l- u) d^, (2.111) 



4 



dt ^:^u{l- u) d^u (1 - u) du. (2.112) 



Expressing the equation of motion in terms of the new time variable u yields 

4 



^u{l-u) duu (1 -u)du + uj^ {q, u) + 4ieu {1 - u) ) x {q, u) = 0. (2.113) 



For convenience we divide by ^ and obtain 

u{l-u) duU (1 -u)du + —co^ [q, u) + ier^u {1 - u) j x (q, u) = 0. (2.114) 

In the next step we use an ansatz for the mode function x {q, u): 

X{q,u)^u"{l-u)%{q,u) (2.115) 
with independent parameters a and . Inserting this ansatz into (2.114) gives after 
^ As a, (3 are constant they can be easily distinguished from the Bogohubov coefficients. 
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a tedious calculation, which is carried out in A.4, the differential equation 

u"+^l-uf^'xd^{duV {q,u)) 
+ (2au°'+^ (1 - - 2/3u"+^ (1 - uf+^ 

(1 - uf-^^ - (1 - uf+^^ X {duV {q, u)) 

+ (a 

+/? - 1) ti°+' (1 - uf - 2a^u'^-^^ (1 - uf^^ 
+Q;ii" (1 - uf-^'' - Pu^^^^ (1 - u)^+^ 
-au''+^ (1 - uf-^^ + (1 - uf^ X 7] {q, u) 

+— {l + ql + ql + ql + 2q^eTU - 2q^eT {I - u) + s^t'^u 

+£^2 (1 - uf - 2e^T\ {l-u)+ 4i£ii (1 - u)) (1 - uf x rj {q, u) (2.116) 

= 0. 

Wc shall simplify this lengthy equation by rewriting it as a hypergeometric dif- 
ferential equation [54]^ 

{u{l - u) dl + {c - {a + b + 1) u) - ab) rj {q, u) = 0. (2.117) 

To this end we make an ansatz for a, b and c 

C— 1 — ITUJo, 

6=1 + -^-'^ + '^. (2.118) 

where we have used 

cu^ = cu^ (m = 0) = 1 + + + (ga - erf , (2.119a) 

a;2 = a;' (li = 1) = 1 + qj + ql + (q^ + erf . (2.119b) 

Using these definitions, the pref actors of 9^77 and i] are given by 

{duV{q,u)) : 

{c-{l + b + a)u) 

= (1 — ITUJo — 2u + ITUJoU — ITUJiU) (2.120) 



^ The hypergeometric differential equation is a special form of Riemann's differential equation. 
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and 

7] {q, u) : 



-ab = -eV + isr^ + ^ + ^ - ^ + ^ _ i^^^. (2.121) 



The hard part is now to check the prefactors in (2.116). We have to divide our whole 
equation by m""*"^ (1 — u)^~^^. Comparing the prefactor yields the coefficients 

« = ^ = — . (2.122) 

Given these expressions, we finally check the equivalence of the terms proportional 
to f] {q, u). The prefactor of 7] {q, u) in (2.116) and the right-hand side of (2.121) have 
to be equivalent. Thus we want to check whether this is the case or not. Writing 
both expression on one side we obtain 



c^o (1 - m) - UJlu^ + l + q^ + 2q^eTU - 2q^eT (1 



u 



I 



^e^tW - 2e^T\ (1 - h) + sV^ (1 - uf + Ae\'^u (1 - u) 

-ujIu-ujIu + ujIu^ ^ujIu^ ^Q. (2.123) 

We pick out all terms with cjq, and find^ 

+ 2ujIu - ulu^ - ulu^ -ul + ul^ -1-^ + 2q^£T - £^t^ - AqzSTU. (2.124) 

Inserting this expression in (2.123) leads to 

e^t'^u'^ - 2e^T\ + 2e^T'^u^ + e^r^ - 2e^T\ + e^t^u^ + Ae^t^u - Ae^t^u^ = (2.125) 

and as all terms vanish we find 

= 0, (2.126) 

which is true and completes the check. 

Our next goal is to derive an expression for the one-particle distribution function 
F {q.,t) in terms of the solutions of the hyper geometric differential equation (2.117). 
This equation has three singular points u = {0, 1, oo}. A convergent solution for 
all ti < 1 can be found for all variables a, b and —c^Nq. As c G C in our case, the 
ansatz 

,^ X „ (a + n)(b + n) 

V{(1:U) = yanU , a^+i = y— — — ra„ (2.127) 
^ (l + n)(c + n) 



— Wq + 2WqU — OJqU^ — OJqU^ — Lof +UJq 

= -1 - ^ + 2q3ST - e^T^ + M + — 2q3STU + e^r'^u -u - - 2q^eTU - e^r'^u 
= -1 - ^ + 2q3ST - e^T^ - 4q3STU. 
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is well defined. For oq = 1 we get the so called hypergeometric function iT\ (a, 6, c; . 
The function 77+ can be identified with the so-called "regular solution" . 

77+(g,M) = 2-Fi(a,6,c;«). (2.129) 

The second linearly independent solution of the hypergeometric differential equa- 
tion is given by 

u^STx (1 + a - c, 1 + 6 - c, 2 - c, -u) = u^'" (1 - uf'"'^ 2^1 (1 - a, 1 - 6, 2 - c, -u) . 

(2.130) 

As a matter of fact, this function can be identified with 

rf {q, u) = -u-^" (1 - u)'"^^ (1 - a, 1 - 6, 2 - c; ?/) . (2.131) 

To get the mode functions (g, u) we have to resubstitute rf^ {q, u) in our ansatz 
(2.115) 

X+ (g, u) = N+u'^ (1 - uf 2^1 (a, c; ^/) , (2.132a) 
X~ (g, xi) = A^-ii"" (1 - u)'^ 2^1 (1 - a, 1 - 6, 2 - c; li) , (2.132b) 

with as normalization constants. These constants are calculated in A. 5 

N+ (q) = =, (2.133a) 

V2a;(g,0) (^ (g, 0) - 713 (g, 0)) ' 

AT- (0) = =. (2.133b) 

V2u;(g,0) (a; (g, 0) + 713 (g, 0)) 

We are now able to calculate the one-particle distribution function according to 
(2.75) 

F{q,t)^2\(3{q,t)f, (2.134) 
as we know each term appearing in the definition of (3 (g, t) 

/3 {q, t) = -ie^K+ {q, t) {dt + icu {q, t)) x+ (g^, t) . (2.135) 

To proceed we transform all functions in (5 (g, t) to the new time variable u. This 

^ It is a special case of the general hypergeometric function pFq{ai, . . . , a^; 61, . . . , 6^; z). Its 
radius of convergence is R = 1 and its serious expansion is given by 



n=0 

where ()„ denotes the Pochhammer symbol. 



2J-i(a,6,c;n) = ^^%^^, (2.128) 
^ (c)„ n! 
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gives 



dt = du {dtu) ^-u{l- u) du, (2.136) 

T 

s]_ (q) = (q, u) - TTg (q, u) = (u (q, u) - TTg (q, u)) (u (q, u) + TTg (q, u)) (2.137) 
and thus^ 



F{q,u) 



uj {q, u) 



iV^ 



X 



u{l — u) du + icu {q, u) ) (1 — u)^ 2^1 (a, c; u) 



. (2.138) 



(++) 

We focus on the term (++) first. Tfie derivative term gives 

du {u°- (1 - uY iT\ {a, b, c; uj^ = au"'~^ (1 - u)'^ 2^1 {a, b, c; li) 

- /3m" (1 - -u)'^"^ 2J^i (a, c; -u) + -u" (1 - — 2 ^"1 (1 + a, 1 + 6, 1 + c; «) , (2.139) 

c 

where we have used the relation 

nh 

(2.140) 

The term (++) can thus be written as 
2 / 

-u (1 — u) \^au"'~^ (1 — u)^ iT\ (a, 6, c; m) — /3m" (1 — uf^^ 2^1 (a, 6, c; m) 
+m" (1 — u)^ — 2-^1 (1 + a, 1 + 6, 1 + c; m) j + iuj (g, m) (1 — u)^ 2^1 («, c; m) 



(x6 

9u2^i (a, i), c; ii) = —2^1 (1 + a, 1 + &, 1 + c; li) . 
c 



(2.141) 



As a, ^ are purely imaginary we get 

m" (1 - uf 



|galn(«)g^ln(l-u) | 



(2.142) 



F{lu)=2\P{q,u)Y 



(f u (1 -u)du + iw) ((f, u) It" (1 - u)'^ 2^1 (a, 6, c; u) 



(1 + ^) |iV+|' (f « (1 -u)du + iw) (1 - uf 2.F1 (a, 6, c; «) 
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The remaining terms will be denoted as 



Fi {q, u) ^-u{l- u) {l + a,l + b,l + c;u), (2.143) 

T C 

F2 {q, u) = {uj {q, u) - {1- u) ujq - uuji) 2 J"i (a, b, c; u) , (2.144) 

such that the one-particle distribution function as a function of the transformed time 
variable u is given by 

F(q,u) = |iV+f (1 + ^^^^) \F,{q,u)+iF2(q,u)\\ (2.145) 

According to the previous discussion of the quasi-particle picture, we are mainly 
interested in the asymptotic behaviour of the one-particle distribution function. 
Hence we investigate F {q. u) in the limit u — )■ 1^ corresponding to t ^ 00. The 
one-particle distribution function for asymptotic times is given by 

F{q,u^l-)^ 

\N+ {q)\' (^1 + (^"'^ ^ 1-) +ii^2 iq,u^ (2.146) 

In A. 6 we explicitly show that this can be simplified to 

2sinh (ttt {2st + luq — Ui) /2) sinh (ttt {2eT — uJo + <^i) /2) 



F{q,u^l-)^ 



sinh (tttcuo) sinh (7rTa;i) 

(2.147) 



2.5 Optimal Control Theory 

In this section we shall give an introduction to optimization followed by a presentation 
of important equations. These relations will be derived step by step, beginning with a 
cost functional and ending with a discussion about the gradient and the optimization 
direction. Finally, the implementation on a computer is discussed. 

Optimal control theory is a mathematical optimization method operating under 
certain control laws. These laws are necessary to have control over the calculation 
and force the system into a state for with an optimality criterion is achieved. To per- 
form these calculations one introduces a cost functional J, which depends on input 
parameters, constraints and control variables. Usually the constraints are inequal- 
ities, chosen such that they force the system away from unwanted configurations. 
The control variables on the other hand are governed by differential equations and 
evaluated such that they do not violate the boundary conditions. The goal is then 
to minimize the cost functional under the given constraints. 

In our case, optimal control theory forces a given potential A (t) into a stationary 
point concerning the asymptotic particle number density. In order to derive the 
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corresponding equations, we start with the differential equation (2.92). Due to the 

fact, that F{q,t), G{q,t) and H {q,t) form a complete system, we have to take 
all of them into account. Furthermore we introduce the Lagrangian parameters 
A*F {q, t) , l^G {q, t) , l^H {q, t) needed for the adjoint differential equation. In addition 
we define the inner products 

(/ (t) , g {t))t = / / (g, t) g {q, t) dq, (2.148) 

Jr 

if, 9)- f {f{t),gmdt. (2.149) 
Jr 

The most important part regarding optimal control theory is to find an appropriate 
cost functional. In our case, we want to maximize the asymptotic particle number, 
so we start with the expression 

J{F{q,t),G{q,t),H{q,t),A)^-^ f F {q,T) dq + p[A], (2.150) 

Jr 

with A(t) as control variable and p[A] being an additional functional. As we try to 
keep the derivation as general as possible we specify the form of p[A] later. 
The next step is to introduce constraints given by 

Xf {q, t)^F {q, t)-W {q, t) G {q, t) = 0, (2.151) 
Xg (g, t)^G (g, t) + W (g, t) F (g, t) + 2u (g, t) H {q, t)-W {q, t) = 0, (2.152) 
A^^ (g-, t)^H (q, t) - 2uj (q, t) G (q, t) = 0, (2.153) 

where as usual 

4(^) = l+??+?2^ (2.154) 

(g-, t)=el + (ga - A, {t)f . (2.155) 
Then we can write down the Lagrange function 

L = J{F, G, H, A) + {Xp, lip) + (Ag, /^g) + (A//, iih). (2.156) 
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We search for a stationary point of the Lagrange function 
L^J + J dqj dt {{updt + Whg) F 

= J + j dq (^hfF \r- J dt {iipF) + j dtWiioF ^^^^ 
+ AtcG \r- J dt {fie) G + J dt {-WiiF - ^cofJ^H) G 
+11hH \r- j dt ifinH) + J dt2uiiGH - j dtW^c 
thus 

which gives the adjoint differential equation 

fiF (g, t) = W (q, t) i^G (g, t) , (2.159) 
fiG (g, t) = -W {q, t) fiF {q, t) - 2uj {q, t) fin {q, t) , (2.160) 
fiH {q, t) = 2cu {q, t) /^g {q, t) • (2.161) 

One has to solve this equation backwards in time with the initial conditions 

^^ = ^,te-r)-Ti0. (2.162) 

^g|^ = MG(«.T)=0, (2.163) 

^^=M«(f.r)i0. (2.164) 

Finally, we proceed with the derivation of the gradient for a given potential. First 
we split the Lagrange function into two parts 

L{F,G,H,iiF,iiG,IJ'H,A)^p[A] + J dq J dtW {hgF - iifG - hg) 

+ 2 [dq [dtcu {-hhG + ugH) + L {F, G, H, fiF, Ac, fin) , (2.165) 



where all irrelevant terms are absorbed into L. 
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Now we define the reduced cost functional from the relevant part 



J (A) = p[A] + JdqJ dtW (i^gF - i^fG - i^g) 

+ 2jdqJ dtcu {-hhG + ugH) (2.166) 

which gives us the gradient 

V J = ^ = 1^ - 2 / cig (-/.^G + ^^GH) ^ 

SA 6A u ^2.167) 

+ / dq-^dtiiiGF - iipG - hg) ■ 

So far we have not specified the functional p[A]. Our idea is to use this term as 
a punishment term to prevent field configurations with higher field strengths from 
dominating the calculations. Therefore, we introduce an upper limit Emax and a 
lower hmit for the field strength Ej^in and choose 

p[A\ = -Ati / dtlii (Emax (t) +A)-fi2 [ dtln (-Emin (t) - A) , (2.168) 

where assesses the importance of the constraints. This fixes the additional 

term for the gradient 

Note that such a term for the constraints still allows us to implement negative field 
strengths. 

During calculation we have to perform a line search to minimize the cost func- 
tional. We proceed by stepwise calculation of the problem 

J (^Ai-aiWJ{Ai)j , (2.170) 

where the gradient — V J {Ai) gives the search direction and ai the step length. A 
crucial point here is to find a^s such that convergence can be achieved. In our case 
we simple use bisection method and apply the so-called Wolfe conditions to the 
algorithm 

J (^Ak - akVJ (Ak)) < J (Ak) - c,ak{VJ (A^) , V J (Ak)), (2.171a) 
(VJ (^Ak - a;feVJ(A)) , VJ(A)) < C2(VJ(Afe) , VJ(Afc)), (2.171b) 
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where < Ci < C2 < 1. The first inequahty imphes that the cost functional becomes 
smaller for each step whereas the second one states that the gradient vanishes at a 
stationary point. 
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3 

Numerics 



Analytical solutions, like the one presented in chapter 2, are only available for a 
very restricted number of field configurations. Hence, one has to rely on numerical 
approaches in general. As we have employed a numerical approach in our calculations 
we have included this chapter for several reasons: It explains the effectiveness of 
the differential equation formalism and introduces several other transformations. 
Moreover, it provides information about the performance of different algorithms, 
shows the results of various numerical benchmark tests and demonstrates the pros 
and cons of discretized potentials. 

3.1 Differential Equation 

As already shown in chapter 2, we benefit from transforming the integro-differential 
equation 



with W (g, t) = el (q) = ml + qf + g|, (g, t) = ml + (q- el(t)) 



for computational reasons. In the following subsections we introduce additional 
transformations and discuss their applicability. 




(3.1) 



into a coupled differential equation of first order 



F{q,t) 
G (q, t) 
H{q,t) 



W{q,t)-G{q,t), 

W (g, t){l-F {q, t)) - 2uj {q, t) ■ H {q, t) , 
2uj {q,t)G{q,t), 



(3.2a) 
(3.2b) 
(3.2c) 
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3.1.1 Transformation 1 



In this subsection we show a transformation, which, we hope, could prevent the 
calculation from producing underflow errors. This should result in a more stable 
code. 

We start again from the intcgro-diffcrcntial equation (3.1). However, we write 
F it) = al {t) this time with a being a simple constant. This gives 



(3.3) 



al {q, t) = W {q, t) / dt'W {q, t'){l- al {q, t')) cos {29 {q, t, t')) . 
Then we introduce the two auxiliary functions 

I dt'W {q,t') {1 -al{q,t')) cos {29 {q,t,t')), 
I dt'W {q,t') {1 -al{q,t')) sin {29 {q,t,t')). 



(3.4) 
(3.5) 



As the system is basically the same as before, the results for the derivatives are 
similar to the previous calculation 



G {q, t)^W {q, t) (1 - al {q, t)) - 2u {q, t) H {q, t) , 
H{q,t) = 2u {q,t)G{q,t). 

We have now an equation for / {q, t) instead of F {q, t) 

I W/a 

= -aW -2u 
\ 2a; 

Using the initial conditions 

I {q,t^ -oo) ^ G{q,t^ -oo) ^ H {q,t ^ -oo) 
we calculate / {q, t) and recover F (g, t) in the end by 

F{q,t) = al {q,t). 



(3.6) 
(3.7) 






fo 


G + 


w 




\0 



(3.8) 



(3.9) 



(3.10) 



The main advantage of this transformation is that / {q, t) is larger than the one- 
particle distribution function for a < 1. 



3.1.2 Transformation 2 

The idea of this transformation is to rewrite the differential equation (3.2) in such 
a way that it is solved on a compact interval. This should help saving computer 
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time and, moreover, the calculations should become more stable as one needs fewer 
evaluation steps. Thus, we start from the coupled differential equation (3.2) and 
employ a transformation, which will be adopted for field configurations of the type 

A{t) = J2 -eintanh (^-^^^ , E{t) = ^e.sech^ (^^^) . (3.11) 

We will distinguish the dominant electric field (index 1) from all other ones (index 
i 7^ 1) until the end of this subsection. The introduced transformation is similar to 
the one used to derive the analytic solution for the single pulse potential in chapter 
2 

t ^ u = ^ (^1 + tanh (J-^^ , (3.12) 

[-00,00] ^ [0,1]. (3.13) 

At that point one gets an impression of the disadvantage of this transformation. The 
most important parameter in the transformation is the pulse length ti, which restricts 
our calculations to setups with a dominant long pulse. The adjective 'dominant' 
means in our notion ti ^ Ti and ei S> Sj. Using the new variable u we obtain the 
following transformed vector potential 

A{u) = (m) + (u) 

i 

^\ ^ /Tiarctanh(2ii- 1) -io,»\ 
= — ri£i (2m — 1) — Ti^jtanh 1^ ^- -j (3.14) 

and transformed electric field 

i 

= AsMl + (1 - t^-h^ (^ ^i^^<^tanh(2n-l)-to, ^| ^^ ^3^^^ 

where io,i is the time lag between the dominant pulse and pulse i. It may look like 
a big disadvantage that the terms in the sum become harder to compute. Analysis 
of the computation time, however, shows that the compression of the time interval 
from [—00, 00] to [0, 1] outweighs the drawbacks. 
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The resulting differential equation looks similar to (3.2) 



duF {q, u) = W {q, u) ■ G {q, u) , (3.16a) 
d^G (g, u) = W {q, u)-{l-F {q, u)) - 2uj {q, u) ■ H {q, u) , (3.16b) 
duH {q,u) = 2uj {q,u) ■ G {q,u) , (3.16c) 

however, the input functions are different 

^ 2ull-t ' ' ^ ylel{q) + {q^-A{u))\ (3.17) 

el{<i)=ml + ql + ql ^i^^^^^^^^' (3-18) 

In the end we want to give some hints concerning numerical efficiency. One can 
observe in (3.14) and (3.15) that the term Tiarctanh {2u — 1) — io,i remains unchanged 
for a given value of u. Therefore it is much faster to calculate this term once and 
read it from the memory instead of calculating it several times. The same reasoning 
applies to r {q, u). 



3.1.3 Low Density Approximation 

Low density indicates that the particle number is low such that F {q, t') <^ 1 for all 
t e R. Accordingly, we approximate the Pauli-blocking term 

(1-F(g,0)~l (3.19) 
such that the integro-differential equation takes the form 

F{q,t)^W{q,t) f dt'W{q,t') cos {29 {q,t,t')). (3.20) 

''Wac 

The benefit of this approximation is that the history of the process itself is ignored, 

which leads to a simple integral. N.B.: This advantage vanishes, however, if one 
deals with the equivalent differential equation 

F{q,t)^W{q,t)-G{q,t), (3.21) 

G (g, t)^W (g, t) - 2uj (g, t) ■ H {I t) , (3.22) 

H{q,t) = 2uj{q,t)G{q,t), (3.23) 

as the neglect of the term {1 — F {q,t')) gives only worse results. 



36 



3.2 Comparison of Transformations 



After deriving three distinct differential equations (2.99), (3.8), (3.16) to describe 
pair creation theoretically, the question of the "best" method arises. Thus we first 
have to identify the favorable transformation and correspondingly a well-adapted 
solver. In our computation we use the configuration 

A{t) = J2 -^^^^ ^^^^ (^V^) ' ^(^) = E (^V^) ' (3-24) 

where we choose specific configurations to investigate three aspects of the different 
transformations: The first and most crucial one is runtime^, the others arc accuracy 
and stability. Usually accuracy is the most important specification, but according 
to the thesis of M. Orthaber [28] we know that a sufficiently good accuracy can be 
achieved. Nevertheless we will also test it on our own, especially with configurations 
that are comparable with the analytical solution. 

In the following we perform runtime and accuracy tests for the different transfor- 
mations and refer to them as 

Standeird differential equation (2.99), 

Stand£ird(I) differential equation (2.99) with an additional implicit solver. 
Transformation 1 differential equation (3.8) with additional parameter a = 10~^, 
Transformation 2 differential equation (3.16), 

Transformation 2(1) differential equation (3.16) with an additional implicit solver. 

We have used a Runge-Kutta 8 solver from Numerical Recipes 3 [55] for non-stiff 
problems and the RADAU5 solver from E. Hairer [56] for stiff differential equations. 
In numerical calculations, we cannot start at asymptotic times t — > — oo nor can we 
end at t — >■ oo. In our simulations, the starting point was chosen such that the vector 
potential was numerically constant up to a relative precision of 10"''', corresponding 
to a vanishing electric field. The end point was reached once the potential becomes 
constant again. 

For the system "Standard(I)" computations using the implicit solver have been 
performed from initial time up to a maximum of 10 per cent of the pulse length of the 
first pulse. In contrast, for the setup "Transformation2(I)" the implicit calculation 
has been done for u starting at 10~^ and ending at 0.1. Then the solver has changed. 

Choosing the right end point for the implicit solver is a very tricky issue. If 
one works too early with an explicit solver, the computation time will increase. If 

^ The values given for the runtime should only hold as an indicator whether a transformation 
works well or not, because for a detailed analysis one would have to perform hundreds of simulation 
runs and analyze all data. 
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one waits too long with a change of the solvers one runs into troubles regarding 
accuracy. Even worse, it could happen that the solver produces completely wrong 
results, although the calculation looks stable ^. 

For the purpose of a systematic evaluation of the various transformations, we 
have picked out a sample of the later used field configurations and compared the 
obtained results. We analyze the runtime performance in Tab. 3.1 and the accuracy 
performance in Tab. 3. 2. Comparing these results one gets the impression that the 
Standard solver is slow, but works for every field configuration under consideration 
as one can see in Tab. 3. 2. The more elaborated solvers "Standard(I)" and "Trans- 
formation2(l)", on the other hand, are the better choice for certain configurations. 
One has to be aware, however, that these solvers cannot be used for every field 
configuration as presented in Tab. 3. 2. 



^ A possible explanation for this effect could be the need of a higher precision for an implicit 
solver. 



38 



Table 3.1: Runtime comparison for different field configurations and 
solvers. The results are normalized to the results of the Standard 
solver to get a better impression of the pros and cons of the various 
transformations. The solver for St (I), Tr2 and Tr2(I) are faster in every 
calculation except those with ten small pulses.^ 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


ot 


Q+/'T^ 
bt(lj 


11-2(1 j 


irz 


TV-1 

in 


£ = 0.1,r = 100^ 


1.000 


0.675 


0.689 


0.910 


2.212 


E — U.UUO, T — O 


i.UUU 


U. / 40 


U.Doz 


U.OiZ 


i.Oii 


£l = 0.1, Ti = 100, 
£2 — U.Ui, T2 — z 


1.000 


0.743 


0.611 


0.893 


2.060 


£1 = 0.1, ri = 100, 

£2 — U.Ui, T2 — -^,10 — >5U 


1.000 


0.674 


0.630 




2.138 


£ = ±0.005,T = 3,io = 53 


1.000 


0.839 


1.706 


1.458 


1.802 


£ = 0.005, r = 3, to = 53 


1.000 


0.827 


1.843 


1.769 


2.203 


£1 = 0.1, Ti = 100, 

Si = ±0.005, Ti = 3, to = 9^ 


1.000 


0.654 


0.422 


0.558 


2.213 


£1 = 0.1, Ti = 100, 

£i = 0.005, Ti = 3,to = 9^ 


1.000 


0.650 


0.430 


0.543 




£1 = 0.1, Ti = 100,to,i = 200, 
£i = 0.005, Ti = 3,^0 = 24.95^ 


1.000 


0.926 


0.306 







^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long 



pulse and ten short pulses 

The tables Tab. 3.1 and Tab. 3. 2 provide much information with regard to coding 
fast and reliable computer program. We can conclude that the use of an implicit 
solver at the beginning of a calculation is a competitive alternative compared to the 
standard solver. The surprising failures are possibly driven by numerical instabilities, 
which could probably be overcome by using a different implicit solver. A good option 
would be a Gauss-Runge-Kutta solver of 4th or higher order. 

These tables further illustrate that the solver "Standard" produces correct results 
for any given configuration above a certain numerical limit, which is set by a combi- 
nation of the field strength parameter and the pulse length. Additionally, we observe 
that the solver "Transformation2" gives the same results as the solver "Standard" for 
all field configurations possessing a dominant pulse, however, performs substantially 
faster in these cases. 
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Table 3.2: Comparison of the particle number for different solvers and 
field configurations. The solver St works for any field configuration, 
whereas the implicit/explicit solver St (I) surprisingly fails for a single 
pulse configuration. The solver Tr2(I), on the other hand, performs 
pretty well for all field configurations except those it was not designed 
for^, i.e. configurations without a dominant pulse. 



Field Parameters 
ie[l], r[l/m], tO[l/m]) 


St 


St(I) 


Tr2(I) 


Tr2 


Trl 


£ = 0.1, r = 100^ 


1.529 


1.668 


1.529 


1.540 


1.529 


e = 0.005, T = 3^ 


5.856 


5.856 


5.856 


5.933 


5.856 


£l = U.i, Ti = iUU, 
£2 = 0.01,r2 =22 


4.Z4U 


4.Z4U 


4.Z4U 


4.Z4U 


A OA n 
4.Z4U 


c — n 1 T- — 1 nn 

£2 = 0.01,r2 = 2,to = 302 


A noi 


zL noi 


A noi 




4 n91 


£ = ±0.005,T = 3,to = 53 


8.700 


8.700 


30.99 


31.03 


8.700 


e = 0.005, T = 3,to = 5^ 


3.782 


3.782 


5.241 


5.267 


3.782 


£l = 0.1, Ti = 100, 

Ei = ±0.005, Ti = 3, to = 9^ 


1.421 


1.421 


1.421 


1.420 


1.421 


£^ = 0.1, Ti = 100, 

Ei = 0.005, Ti = 3,^0 = 9^ 


1.353 


1.353 


1.353 


1.353 




£^ = O.Lri = 100,to,i = 200, 
Ei = 0.005, = 3,to = 24.95^ 


2.259 


2.259 


2.259 







^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long 



pulse and ten short pulses 

In order to conduct stability tests, we have performed calculations upon changing 
the parameters only slightly. In this way we can identify parameter regions in which 
our numerics work properly. These calculations indicate that the particle distribution 
function is more difficult to calculate for negative values of 53 than for positive ones. 
Probably this happens, because particle pair creation occurs for negative values of 
53 at earlier times such that one woTiId have to start the calculations at earlier 
initial times, too. Another aspect worth mentioning is that a change t — )■ — t in the 
potential gives a mirrored distribution function F {q^jt 00) F {—q^^t 00). 
This is a further monitoring option for our solvers and in perfect agreement with 
ideas presented in C. K. Dumlu et al. [44]. 
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3.3 Full Momentum Space 



In the previous section we only investigated the performance of the various differential 
equations. In this section we shall give ideas for further optimization, especially for 
the calculation of the particle spectrum in full momentum space. An option, for 
example, is to use an interpolation algorithm avoiding calculations of less important 
data points. Another idea would be to use an improved integration routine such as 
Romberg integration. 

First we investigate the issue of interpolation routine. To test them we have used 
single pulse configurations involving configurations for qi — q2 — 0. For the purpose 
of testing we use the differential equation (3.16), but the results remain valid for the 
other transformations as well. 

We tried a ID cubic spline interpolation [55] in direction of fixed q^. On the other 
hand, we use a linear spline interpolation [55] in the direction of fixed q± as well as in 
case of errors in direction of fixed q^. As a matter of fact, we have also tested other 
routines, but they all suffer from severe disadvantages. A 2D cubic spline for example 
does not work properly, because of the typically irregular data grid. Working with 
such a method would cause a lot of additional computer time without the certainty 
that the result would be worth it. A 2D linear spline interpolation as well as Laplace 
interpolation are not able to produce accurate results as higher order interpolation 
routines are a necessity at least in direction of fixed q^ . 

We investigate the performance for calculating the particle number for various 
parameters. It can be seen from the data in Tab.3.3B that the results for Aq± G 
[0.1, 0.5] m and Aq± = Aq^ are in good agreement with the analytical calculation. As 
expected, the results get better the more data points are involved. Tab.3.3A, on the 
other hand, shows that the runtime increases significantly by decreasing Aq^ — 0.05m 
to Ag3 = 0.01m whereas the accuracy of the results is not improved. 

A similar analysis is performed for a different parameter set in Tab. 3. 4. These 
results indicate that we have to use a step size of at least Aq± = 0.1m for these 
parameters. A decrease in the momentum spacing is, however, accompanied by 
an increase in runtime. As a matter of fact, if we choose the spacing too small, 
the runtime increases substantially whereas the accuracy is not affected anymore. 
Accordingly, finding the appropriate value for Aq is an important but nontrivial 
task. 
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Table 3.3: Runtime comparison for different momentum spacings 
for a single pulse in particle distribution with e = 0.005, r = 3m~^ 
in A. The results are normalized to the fastest calculation and 
all momenta are given in units of the electron mass q[m]. In B 
the results for the particle number are normalized to a calculation 
based on the analytical result. 



A 


Aq^ = 0.5 


Aq^ e [0.1,0.5] 


Aq± = Aga 


Important q^^ 


1.000 


1.631 


3.063 


Aqs = 0.05 


1.210 


1.621 


3.263 


Aqs = 0.01 


3.810 


6.968 


62.705 


^ Aqs = 0.05, 10% 
cally 


Tolerance, whether a Aq^ is calculated 


is chosen dynami- 


B 


Aq^ = 0.5 


Aq^ e [0.1,0.5] 


Ag_L = Aga 


Important q^^ 


1.0255 


1.0021 


0.9979 


Aqs = 0.05 


1.0253 


1.0039 


0.9997 


Aga = 0.01 


1.0245 


1.0047 


1.0000 



^ Aqs = 0.05, 10% Tolerance, whether a Aqs is calculated is chosen dynami- 
cally 
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Table 3.4: Runtime comparison for different momentum spacings 
for a single pulse with e = 0.1, r = 100m~^ in A. The results for 
the particle number are normalized to the fastest calculation and 
all momenta are given in units of the electron mass g[m]. In B 
the results for the particle number are normalized to a calculation 
based on the analytical result. 



A 
A 


Aq^ = 0.5 


Ag^ G [0.1,0.5] 


/\q± — Aga 


imporxanx 


1.000 


2.129 


0.4101 


Aqs = 0.05 


1.345 


2.943 


4.843 


Aqs = 0.01 


6.235 


14.058 


98.640 


^ Aqs = 0.05, 10% 
cally 


Tolerance, whether a Aq^ is calculated 


is chosen dynami- 


B 


Ag^ = 0.5 


Ag^ e [0.1,0.5] 


Ag_L = Aga 


Important q^^ 


0.771168 


1.02592 


1.02585 


Aq3 = 0.05 


0.770376 


1.0005 


1.0005 


Aqs = 0.01 


0.74475 


0.9998 


1.00086 



^ Aqs = 0.05, 10% Tolerance, whether a Aqs is calculated is chosen dynami- 



cally 

We now turn our attention towards integration algorithms. For configurations 
with g_|_ = it was sufficient to use the trapezoidal rule for integration, but for 
calculations in two dimensions this is not true any more. We have used Romberg in- 
tegration instead to obtain better results. The problem with this integration routine 
is, however, that we need data points outside of the calculated grid. In order to avoid 
calculation of extra data points, which would increase the runtime tremendously, we 
have used a 2D polynomial interpolation [55] instead. As one can see in Tab.3.5A, 
the trapezoidal rule is still competitive for short pulses. For the long pulse, however, 
we need a much smaller step size of Aga = 0.01m in the trapezoidal rule integration 
compared to the Romberg integration. 
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Table 3.5: Comparison of the accuracy when adopting different integra- 
tion routines for a single pulse with A) £ = 0.005, r = 3m~^ and B) 
£ = 0.1, T — 100m~^. In addition the influence of different momentum 
spacings has been investigated. The results for the particle number are 
normalized to a calculation based on the analytical result. All momenta 
are given in units of the electron mass q[m\. 



A 

A 


Aq± 


= 0.5 


Aq± e [0.1,0.5] 


Aqi_ 


= Aq3 




TR^ 


RI^ 


TR 


RI 


TR 


RI 


Important q^^ 


1.0255 


1.0210 


1.0021 


1.0070 


0.9979 


1.0027 


ZA53 — u.uo 


1.0260 


1.0253 


1.0110 


1.0039 


1.0068 


0.9997 


Aqs = 0.01 


1.0226 


1.0245 


1.0038 


1.0047 


0.9997 


1.0000 


^ Agg = 0.05, 10% 
^ Trapezoidal rule ^ 


Tolerance, whether a Aq^ is 
Romberg integration 


calculated is 


chosen dynamically 


B 


Aqi_ 


= 0.5 


Aq^ e [0.1,0.5] 


Aqs_ 






TR^ 


RI^ 


TR 


RI 


TR 


RI 


Important q^^ 


1.2967 


0.7711 


1.0259 


1.0259 


1.0258 


1.0258 


Aq3 = 0.05 


0.7703 


0.7703 


1.0263 


1.0005 


1.0262 


1.0005 


Aqs = 0.01 


0.7458 


0.7447 


1.0010 


0.9998 


1.0008 


1.0008 



^ Aqs = 0.05, 10% Tolerance, whether a Aq^ is calculated is chosen dynamically 
^ Trapezoidal rule ^ Romberg integration 
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3.4 Comparison with Lookup- Table 



We finally want to check whether simulations with precalculated values of the vector 
potential and the electric field yield the same results as calculations where these 
values are calculated on-the-fly. This investigation is motivated by the request to 
calculate the pair production process for arbitrary field configurations as well. 

Therefore, we recompute all field configurations, which have previously been dis- 
cussed in Tab. 3.1 and Tab. 3. 2. In comparison, we have now precalculated the vector 
potential and the electric field for a different number of sampling points. The values 
for a required time step have been calculated by interpolation routines. 

Table 3.6: Runtime comparison for different sets of sampling points. Po- 
tential and electric field have been discretized first. In order to calculate 
intermediate values linear interpolation is used. The results for the runtime 
are normahzed to the results of the "Standard" solver. 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


St 


Trans2(I) 


512 
Points 


4096 
Points 


32768 
Points 


£ = 0.1, r = 100^ 


1.000 


0.689 




0.784 


2.532 


£ = 0.005, r = 3^ 


1.000 


0.682 


0.331 


4.470 


22.82 


£^ = O.Ln = 100, 
£2 = 0.01, r2 = 22 


1.000 


0.611 


0.202 


0.527 


1.900 


£i = 0.1, ri = 100, 

£2 = 0.01, r2 = 2,to = 302 


1.000 


0.630 


0.220 


0.526 


1.872 


£ = ±0.005,r = 3,to = 53 


1.000 




0.099 


1.406 


4.441 


£ = 0.005,T = 3,io = 53 


1.000 




0.051 


1.535 


5.016 


£i = 0.1, n = 100, 

Ei = ±0.005, n = 3, to = 9^ 


1.000 


0.422 




0.167 


0.699 


£i = 0.1, Ti = 100, 

€i = 0.005, Ti = 3,to = 9^ 


1.000 


0.430 


0.090 


0.143 


0.692 



^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long pulse 
and ten short pulses 
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Tabic 3.7: Comparison of the particle number for different sets of sampling 
points. In order to calculate intermediate values linear interpolation is used. 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


St 


Tr2(I) 


512 
Points 


4096 
Points 


32768 
Points 


£ = 0.1, r = 100^ 


1.529 


1.529 




1.586 


1.529 


£ = 0.005, r = 3^ 


5.856 


5.856 


5.733 


5.859 


5.857 


£^ = 0.1, n = 100, 

£2 = 0.01,r2=22 


4.240 


4.240 


2.218 


2.748 


4.212 


£-^ = 0.1, n = 100, 

£2 = 0.01, r2 = 2,to = 302 


4.021 


4.021 


1.672 


2.572 


3.993 


e = ±0.005, r = 3, to = 5^ 


8.700 




8.013 


8.697 


8.701 


e = 0.005, r = 3,to = 5^ 


3.782 




3.435 


3.780 


3.783 


= 0.1, Ti = 100, 
Si = ±0.005, Ti = 3,to = 9^ 


1.421 


1.421 




1.006 


1.413 


£^ = 0.1, n = 100, 

£i = 0.005, Ti = 3,to = 9^ 


1.353 


1.353 


7.925 


9.193 


1.345 



Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long pulse and 



ten short pulses 

The number of sampling points is very important as one can see in Tab. 3. 6 and 
Tab. 3. 7. Interestingly, we can observe a lack of accuracy in case of 512 sampling 
points. For most configurations, at least 4096 samphng points are essential. For 
very difficult and complex pulse structures even 32768 points produce an error of 
1%. Additionally, these tables demonstrate that there are field configurations with a 
comparatively long runtime. This is quite surprising as lookup tables arc usually used 
to accelerate calculations. It seems that error correction during a calculation slows 
down the whole computation. This also indicates that there are problems concerning 
stability when using discretized potentials in conjunction with linear interpolation. 

In order to remedy the difficulties with linear interpolation we switch to a more 
advanced interpolation routine, namely cubic interpolation. This interpolation rou- 
tine is more advantageous when dealing with a high number of sampling points as one 
can see in Tab. 3. 8 and Tab. 3. 9. For 512 sampling points there are empty cells due to 
a very long runtime, however, the larger the set of sampling point get the better the 
calculations become. For 32768 sampling points, for example, cubic interpolation 
outperforms linear interpolation as presented in Tab. 3. 8. 
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Tabic 3.8: Runtime comparison for different sets of sampling points. In 
order to calculate intermediate values cubic interpolation has been used. 
The results for the runtime are normalized to the results of the "Standard" 
solver. 



Field Parameters 
[e[i\, T[i/mJ, tU[i/mJj 


St 


Trans2(I) 


512 
Points 


4096 
Points 


32768 
Points 


— U. i . / — iUU 


i .()(}(} 


U.boy 


u.yt) io 


u.yoo 


1 111 
i . i i i 


£ = 0.005,T = 3^ 


1.000 


0.682 


1.183 


1.585 


5.734 


El = 0.1, Ti = 100, 
£2 = 0.01,T2=22 


1.000 


0.611 




0.7291 


0.9519 


£^ = 0.1, n = 100, 

C2 = 0.01. 72 = 2. to = 30- 


1.000 


0.630 




0.8945 


1.167 


e = ±0.005, T = 3,^0 = 5=^ 


1.000 




0.2009 


0.2319 


0.4572 


£ = 0.005,T = 3,to = 5^ 


1.000 




0.1824 


0.2396 


0.5574 


= 0.1, Ti = 100, 
Ei = ±0.005, Ti = 3, to = 9^ 


1.000 


0.422 




0.3086 


0.3273 


El = 0.1, Ti = 100, 

Ei = 0.005, Ti = 3,^0 = 9^ 


1.000 


0.430 




0.3076 


0.3311 



^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long pulse 
and ten short pulses 
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Table 3.9: Comparison of the particle number density for different sets of 
sampling points. In order to calculate intermediate values cubic interpolation 
has been used. 



Field Parameters 
(e\l\ tH/ui] tOU/ui]) 


St 


Tr2(I) 


512 
Points 


4096 
Points 


32768 
Points 


£ = 01 r = 100^ 

\j • J- ^ 1 j-yjyj 


1.529 


1.529 


1.45 


1.519 


1.528 


e = 0.005. - = 3^ 


5.856 


5.856 


5.617 


5.838 


5.85 


El = 0.1, Ti = 100, 
£2 = 0.01,r2 = 22 


4.240 


4.240 




3.162 


4.221 


£l = 0.1, Ti = 100, 

£2 = 0.01,T2 = 2,io = 302 


4.021 


4.021 




2.970 


4.00 


e = ±0.005. - = 3. Jo = 5'^ 


8.700 




7.990 


8.631 


8.692 


£ = 0.005,T = 3,io = 53 


3.782 




3.637 


3.791 


3.784 


£l = 0.1, Ti = 100, 

Ei = ±0.005, n = 3, to = 9^ 


1.421 


1.421 




1.127 


1.415 


El = 0.1, Ti = 100, 

Ei = 0.005, Ti = 3, to = 9^ 


1.353 


1.353 




1.04 


1.347 



^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long pulse and 
ten short pulses 
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4 

Results 



In this chapter we present our new findings regarding particle pair production, with 
a focus on the following key aspects: First an analysis of combinations of two pulses 
with two different time scales is done. These combination of fields have been dis- 
cussed recently, but we provide an accurate perspective on setups with new combi- 
nations and also on the full momentum spectrum in comparison to these previous 
investigations [28,32,34]. Secondly an investigation of many-pulse configurations is 
done. In addition, the differences between fields consisting of equally and alternating 
signed pulses are examined. These investigations should produce novel insights for 
an either direct or indirect measurement of the Schwinger effect. Adopting optimal 
control theory, we finally present results on pulse shaping in particle pair production 
for the first time. 
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4.1 Dynamically Assisted Schwinger Effect 



As mentioned previously, wc focus on combinations of pulses with different time 
scales. To this end, we employ the electric field: 

E{t) = J2 ^.sech^ (^^) • (4-1) 

Before ultimately coming to the results, we shall introduce the Keldysh parameter 
7. It is defined for a single pulse as 

7=^ = ^ = —, (4.2) 

where r is the pulse length and tt is the characteristic time for a tunneling process 
from the Dirac sea to the continuum. The Keldysh parameter separates two dif- 
ferent regimes. For a parameter value of 7 <^ 1 the pair production process shows 
characteristics of a tunneling process and is therefore associated with the Schwinger 
effect. On the other hand, a parameter value of 7 ^ 1 indicates that a high number 
of photons is absorbed. This parameter region is referred to as multiphoton regime 
as the multiphoton absorption process dominates the pair production in this case. 
It is the idea of the dynamically assisted Schwinger effect to combine pulses in the 
different parameter regimes. In order to describe such a process properly, we intro- 
duce the combined Keldysh parameter 7c, which is defined based on properties of 
both pulses [32,34]. 

ic^^—. (4.3) 

meT2£l 

Here r2 is the pulse length of the short pulse and ei the maximal field strength of the 
long pulse. This enables us to characterize the field pulses in terms of one parameter 
as illustrated in the following figures. 

In Fig. 4. 1 we investigate the enhancement of particle production when combining 
two pulses, the first one with 71 <^ 1 and the second one with 72 » 1. We further 
define enhancement as the ratio 

Nr 

where Ni, N2 denote the particle number in a single pulse of given parameters, 
whereas Nq denotes the particle number for combined electric fields. As can be seen 
in Fig.4.1 the enhancement becomes largest at a combined Keldysh parameter of 
7c; f=s 2. At this point the pair production rate of the two single pulses is of the same 
order. Remarkably, there seems to be no qualitative difference between calculations 
based on the full momentum space Fig.4.1a-d and calculations based only on the 
parallel momentum Fig.4.1e. 
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1 2 3 4 1 2 




Combined Keldysh Paiameter 



Figure 4.1: Enhancement of the particle number density due to dynamical assis- 
tance. The calculations in a — d are based on full momentum space calculations 
whereas the calculation in e is done in the parallel direction only. The offset be- 
tween the two pulses is to = 0. The absolute field strength is varied in each plot 
line. The same pulse length ri is used in a, c and 6, respectively, whereas the 
field strength ratio 62/61 is the same in a, h and c, respectively. 
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12 3 4 



Combined Keldysh Parameter 

Figure 4.2: Enhancement of the particle number density based on a full momentum 
space calculation. A short pulse is assisting a long pulse with an offset of to = 0. 
The field parameters of the short pulse are varied while the parameters for the long 
pulse are fixed at ei = 0.1, ri = 102. 2m~^. 

In the next calculation, illustrated in Fig. 4. 2, we varied the field strength of the 
short pulse, while keeping the parameters for the long pulse fixed at ei = 0.1, ti = 
102. 2m~^. Again, the effect from the dynamical assistance of photon absorption is 
most pronounced at 7c ~ 2. The highest relative enhancement effect can be seen 
for the configuration 62 = 0.005. On the other hand, the highest absolute values are 
obtained for e = 0.025. The increase in the absolute value of the enhancement for 
increasing values of £2 is simply explained with the higher overall field strength in 
combining the two pulses. 
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rl=50[l/m] 






0.0 0.5 1.0 1.5 2.0 2.5 

Combined Keldysh Parameter 

Figure 4.3: Enhancement of the particle number density based on a calculation in 
the parallel direction only. The offset between the two pulses is to = 0. The field 
strengths are ei = 0.1 for the long pulse and 62 = 0.01 for the short pulse. 

Subsequently, we kept the field strength constant at ei = 0.1, 62 = 0.01 and 
varied just the pulse lengths. The results obtained from these calculations are shown 
in Fig.4.3. Interestingly, all lines meet at 7c > 2.2 whereas the influence of the 
different pulse lengths can be observed in the region < 7c < 2. The limit 7c — >■ 
corresponds to a situation where the first pulse is in the regime 7(7 <^ 1 whereas the 
second pulse becomes adiabatically constant. Thus the outcome corresponds more 
or less to the outcome obtained from a single pulse with pulse strength ei +62- In 
the limit 7c > 2.2 the particle production due to the second pulse 72 ^ 1 starts 
to dominate. In this limit, the effect of the first pulse 71 ^ 1 becomes effectively 
negligible. 
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Figure 4.4: "©MMaMttli^lii the particle distribution fuiffff!6tfTS!ttg?Ja%y combining 
two pulses with different pulse scales. We assume ei = 0.25, ri = 50m^^ for the 
long pulse. The field strength of the short pulse is 62 = 0.025 whereas different time 
scales T2 = lm~^ or T2 = 5m~^ have been used. The offset between the two pulses 
is to = 0. The dashed line is the simple sum of the long and the shorter second 
pulse. In a the two fields strengths have the same sign whereas they have opposite 
sign in b. 



Finally, we also investigate the particle distribution in Fig.4.4. We have used 
a long pulse with parameters ei = 0.25, ti = 50m~^ combined with a short pulse 
62 = 0.025, T2 = lm~^ and 62 = 0.025, T2 = 5m~^, respectively. For a single pulsed 
field 71 <^ 1 no oscillations in the distribution function can be observed. This is 
even true for configurations of two pulses when T2 is not too short. For the even 
shorter multiphoton pulse T2 = lm~^, however, the distribution is deformed leading 
to more than one local maximum. Moreover an interesting behaviour of the particle 
distribution is observed when changing the sign of the second pulse 62 — )■ —£2- This 
effect will be discussed in more detail in the following sections. 
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So far we have seen that one additional short pulse can significantly change the 
distribution and the number density of the created particles. We proceed by in- 
vestigating configurations with more than one short pulse. To be specific, we will 
investigate dynamical assistance of a long pulse by a train of short pulses as well as 
the resulting interference pattern in the momentum distribution. 

We present results on the momentum distribution for a configuration with one 
long pulse £1 = 0.1, ti = 150m~^ and ten identical small pulses with 62 = 0.05, T2 = 
5m~^ separated by a time lag of 10m~^ in Fig. 4. 5. Note that this is only a shorthand 
notation. Explicitly, it means that the first small pulse has a time lag of to,i = — 4.5to 
and the last small pulse has a time lag of to, 10 = 4.5to- Moreover, we investigate 
configurations where the small pulses have the same sign as the long pulse (denoted 
as parallel) and alternating signs (denoted as alternating), respectively. In fact one 
can observe the effect due to each small pulse in the parallel configuration in Fig. 4. 5a. 
For the alternating configuration, however, only the pulses parallel to the long pulse 
yield an enhancement as seen in Fig.4.5b. Again, this effect is due to the dynamical 
assistance, as each photon absorption boosts the particle production only in a small 
momentum range. Furthermore, we observe that the two different configurations 
differ by the range occupied in momentum space. This might be due to the fact 
that all small pulses act accelerating in the parallel configuration whereas half of 
the pulses also act decelerating in the alternating configuration. Finally, the particle 
number at its local minima shows different characteristics as well: They increase in 
the center for the parallel configuration whereas they decrease in the center for the 
alternating configuration. Unfortunately, we could not find a simple explanation for 
this effect. 
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T-i J I- TI Pvallel Mopientum[m] „ , . Paralle.1 Momeiitum[m] 

r IGURE 4.5: Jr^article distribution function for a connguration with a long pulse being 
assisted by ten evenly spaced short pulses. Their field parameters are ei = 0.1, ri = 
150m~^, \e2\ = 0.05, T2 = 5m~^ and the time lag is chosen to be to = 10m~^. We 
present the parallel configuration in a and the alternating configuration in b. 



55 



In Fig.4.6 we present the results for a setup with one long pulse ei = 0.1, ti = 
100m~^ and ten short pulses 62 = 0.005, T2 = 3m~^. Again, the ten short pulses 
are evenly spaced, however, the pulse train is only switched on when the long pulse 
reaches its peak value. Thus, to,i = and to,j > is true. Consequently interferences 
can be seen in the particle distribution function also for higher parallel momenta in 
Fig.4.6. 




02468 10 02468 



Parallel Momentum[m] Parallel Momentum[m] 

Figure 4.6: Particle distribution function for a configuration with a long pulse 
being assisted by ten evenly space short pulses. The characteristic field parameters 
for the long pulse are ei = 0.1, ri = 100m~^ and for the 10 evenly spaced short 
pulses they are £2 = 0.005, T2 = 3m~^. The first small pulse has a time lag of 
to,i = 0. Again we present the parallel configuration in a and the alternating 
configuration in b. The results are compared with a field consisting of all pulses 
with positive field strengths placed on top of each other as shown by the gray line. 



4.2 Interferences in Many Pulse Configurations 

In this section we want to analyze the ideas presented in [43] within the quantum 
kinetic approach. In this paper, it has been suggested that setups with alternat- 
ing fields show characteristics of a Ramsey interferometer. We have performed the 
corresponding computations and additionally examined configurations consisting of 
equally signed pulses. 

The particle distribution function for both cases is shown in Fig. 4. 7. We consider 
= 10 equally spaced pulses with parameters e = 0.02, r = 3m^^ and a specific 
time lag. Most notably, the two configurations show a very different behaviour: In 
the parallel configuration we obtain a width in momentum space corresponding to A^^ 
times the particle distribution function of a single pulse with the same parameters 
as shown in Fig. 4. 7a. For the alternating configuration, however, one obtains the 
corresponding maximum magnitude as shown in Fig.4.7b. Note, that the exact form 
of the distribution function strongly depends on the choice of the time lag to- These 
findings further support the idea of choosing pulses with an adjusted time lag to 
enhance the pair production probability substantially. 
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Similarly to [43] we have also calculated the peak value of the particle distribution 
function for an alternating field configuration consisting of = 2 short pulses with 
parameters e = 0.1, r = 25m~^. The results are illustrated in Fig.4.8 and are 
consistent with the previous discussion. 
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Figure 4.7: Positive interference effects in particle distribution for = 10 short 
pulses with field parameters are e = 0.02, r = 3m~^. For comparison, we plot N'^ 
times the particle distribution function of a single pulse with the same parameters. 
We present the parallel configuration in a and the alternating configuration in b. 
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Figure 4.8: Behaviour of the peak value of the particle distribution function for 
N = 2 short pulses with parameters ei = 0.1, ti = 25m^^ and £2 = —0.1, T2 = 
25m~^. The maximal value of the distribution function is shown as a function of the 
time lag and compared with iV^ times the peak value of the distribution function 
of a single pulse with the same parameters. 
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4.3 Oscillations in Particle Density 



Based on the discussion of the previous section, it has been our main goal to combine 
interference effects with the dynamical assistance effect. Obviously, it is unreahstic 
to do calculations for the whole configuration space, so we have computed only the 
most intuitive setups. Accordingly, we will examine the dynamical assistance of 
a long pulse with a train of short pulses in both the parallel and the alternating 
configuration. Finally, an additional antisymmetric configuration has been used for 
comparison. This setup has been tweaked such that one has only parallel pulses in 
the beginning and anti-parallel pulses in the end. 

The results presented in Fig. 4. 9 exhibit several surprising details of electron 
positron pair production. In this investigation, we defined a maximal field strength 
£max = 0.1 and a maximal energy of the configuration En = 1.33m to avoid that field 
strength effects spoil the result. Moreover, we consider pulse trains consisting of 10 
equally spaced short pulses. The thick blue line shows the particle number for the 
superposition of a single long pulse with the pulse train whereas the red line shows 
the particle number for the pulse train only. 

Wc first focus on the results of the pulse train without a long pulse. We observe 
that the particle number oscillates as function of the time lag. Most notably, big 
local maxima and comparatively less pronounced local extrema arise, at least for 
the parallel and the alternating configuration. In the anti-symmetric configuration, 
however, a double peak structure appears and the local extrema become much more 
pronounced. It is quite remarkable, however, that the oscillatory behaviour as a 
function of the time lag is very regular and stable in all possible configurations, at 
least for large values of to- 

Adding a single long pulse to the pulse train, we find in Fig.4.9 that the pattern 
of local extrema due to the pulse train is inherited by the full solution. We observe, 
however, that a long pulse can keep the small pulses from causing oscillations in 
the distribution function for rather large values of T2. For decreasing values of T2, 
however, the oscillations are damped for higher time lags but stay observable. It 
is also remarkable that the superposition of a long pulse results in a delay of the 
oscillation frequency. Thus, the surprising effect arises that the additional long pulse 
can both increase and decrease the particle number depending on the time lag to- 
Based on Fig.4.9, we may observe interesting aspects. We start our analysis with 
a detailed description of the results for the time interval to e [0, 3m~^] and have a 
look at the pattern of local minima in Fig.4.10. These results indicate that there is 
a strong correspondence between the number of the local minima and the number of 
pulses used. Besides comparing the behaviour of the particle distribution function 
for the first two time intervals in Fig. 4. 11, a more detailed comparison between 
the parallel and the alternating configuration is shown in Fig.4.12. This surprising 
result raises the question whether there is a deeper connection between these two 
configurations. In Fig.4.13, we finally resolve that the local maxima pattern in a 
pulse train can cause oscillations in setups with an additional long pulse. 
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i'lGURE 4.9: rarticle number as mnction oi the time lag Iq between A* = lU 
short pulses with pulse length T2 assisting a single long pulse. The thick blue line 
shows the particle number for the superposition of a single long pulse with the pulse 
train whereas the red line shows the particle number for the pulse train only. The 
maximal field strength is limited to Emax = 0.1 and the maximal energy to 1.33m. 
We present results for the parallel configuration in a, b and for alternating fields 
in c, d. In e, / we started with positive fields and ended with negative fields. For 
a better orientation the values of 1, 10, 100 times the particle number of a single 
small pulse as well as the results for asymptotic time lags are shown. 
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1 2 3 4 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 

Time Inlerval[l/m] Time Inlerval[l/m] 

Figure 4.10: Pattern of the local extrema in the particle number for a different 
number of pulses in the pulse train. A different number of pulses produces a different 
number of local maxima between the peak values. We show the results for pulse trains 
with parameters e = 0.005 and r = 3m~^ in the parallel configuration. 




Time Interval[a.u .] 



Figure 4.11: Damping of the local extrema pattern in the particle number upon 
changing the time interval. 'First interval' refers to a time lag of to ^ [0,3] whereas 
'second interval' refers to a time lag of to G [3,6]. We show the results for a pulse 
train of iV = 10 pulses with parameters e = 0.005 and r = 3m~^. The height of 
the oscillations decreases for higher time intervals. 
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e= ±0.005 
e=+0.005 




1 2 3 

Time Interval[l/m ] 

Figure 4.12: Dependence of the particle number on the used field configuration. 
Again, we show the results for a pulse train of = 10 pulses with parameters 
l^l = 0.005 and r = 3m~^. The red line corresponds to a parallel field configuration, 
while the blue line corresponds to the alternating configuration. 
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Figure 4.13: Inheritance of local oscillations in the particle number from a pulse 
train of = 10 pulses with parameters e = 0.005 and r = 3m~^ compared with a 
single long pulse in addition to the same pulse train. 



61 



4-3.1 Time-scale Comparison 

Similar to the last section we focus on the influence of the chosen configuration on the 
particle number and provide further data for pulse trains in parallel and alternating 
configurations. Again, we use pulse trains consisting of = 10 equally space pulses 
and vary their pulse length T2. 

It can be observed in Fig.4.14, that the oscillations in the parallel case become 
smaller in magnitude the longer the pulse gets. For the alternating configuration, 
however, additional peaks arise which become even more pronounced for longer 
pulses. For long pulses it even looks as if the frequency of the oscillations has doubled 
in the alternating case. On the contrary, the weak local extrema in the symmetric 
configuration seem to vanish. 
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Figure 4.15: Behaviour oi local fluctuations m the particle number with increasing 



pulse length for parallel configuration. The used field strength has been e 
The pattern of the local extrema vanishes for longer pulse lengths. 



0.02. 



We further investigate the parallel configuration and focus on a time lag of about 
to ~ 20. The first interesting observation, as illustrated in Fig. 4. 15, is that the 
small local extrema vanish entirely for longer pulses. The main maximum, however, 
becomes broader. 
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An additional local extremum arises for longer pulse lengths. 
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In Fig.4.16 we have done the same detailed analysis of local extrema, but this 
time for the alternating configuration. The main difference compared to the parallel 
configuration is that the local extrema do not vanish for longer pulse lengths. Instead 
another local extremum raises for longer pulses. Therefore, one gets the impression 
that the rate of oscillations has doubled. This gives rise to a distinctive difference 
between the parallel and the alternating configuration. 
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4-3.2 Field Strength Comparison 

In the last part of this section, we investigate the impact of the field strength on the 
particle number for the parallel configuration. In this case, we use a pulse train with 
= 10 single pulses with pulse length r = 3m^^ and different field strengths 62- It 
can be seen in Fig. 4. 17 that there is a qualitative change in the particle number for 
increasing field strengths. 

Again, we further investigate the details in Fig.4.18. We can see that the shape 
of the particle distribution function shows significant distinctions for different field 
strengths. 
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Figure 4.17: Particle number for a puf 
T = 3m~^ and increasing pulse strength. 
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4.4 Pulse Shaped Fields 



In this section, we explain the basic aspects of pulse shaping. The goal of these cal- 
culations is to maximize the particle number by deforming the electric field according 
to the optimal control algorithm. At the beginning, we introduce the procedure and 
then go ahead and present the data. First, we define an initial electric field, which 
has the form 

E{t) = e sech^ (t/r) . (4.5) 

Then we discrctizc the vector potential and the electric field with 4096 non-equidistant 
points. These sampling points are obtained from a set of equidistant sampling points 
by means of an atanh mapping. This transformation needs to be performed in order 
to obtain the necessary accuracy later on. Additionally, we will compute the gradient 
for the search direction at these sampling points as well. Subsequently, we solve both 
the standard differential equations of first order (2.99) and the adjoint differential 
equation (2.161). We do not only save the results at asymptotic times, but keep the 
values at all sampling points. The computation of the gradient is straightforward 
(2.167) and gives us the search direction via steepest descent method 

Ai+i{t) = Ai{t) - aiVJ {Ai{t)) ■ sech' {t/r) . (4.6) 

The sech^ (t/r)-term is introduced to ensure sufficient damping of the gradient so 
that the electric field vanishes at asymptotic times for every iteration step. Another 
crucial task is to fix the constraints for the optimal control algorithm. Generally 
spoken, the electric field strength increases ad infinitum if there is no upper bound. 
To ensure that this will not happen in our simulations we introduced bounds for the 
field strength. 

The easiest way to do this is to define a maxinmm field strength for all times. 
However, we also tried to use time-dependent constraints. In Fig. 4. 19 we show the 
deformation of the initial pulse s = 0.01, r = 5m~^. In this calculation the upper 
hmit is set at £ = 0.01, r = 6m~^ whereas the lower hmit is set at £ = 0.009, r = 
4m~^. Both hmits hold with a small offset to ensure that there is free space to 
explore. 

We investigated the impact of weaker constraints, too. The corresponding results 
are shown in Fig.4.20. It can be seen that the preferred substructure appearing after 
optimization has periodic form. 

In order to understand this pulse form we have further investigated setups with 
two different initial conditions in Fig.4.21. These initial parameters are e = 0.01, r = 
3m~^ in Fig.4.21a and e = 0.1, r = 20m^^ in Fig.4.21c. In both cases the upper 
limit for the field strength is specified at 101 percent of the initial field strength. As 
one can see from Fig.4.21b and Fig.4.21d, the substructures favored by the optimal 
control algorithm are very similar. Most notably, both the oscillation rate and the 
peak height of the substructure resemble each other. Accordingly, one may conclude 
the existence of an optimal field configuration. 
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According to Fig.4.21 we conclude that the optimal setup shows substructures of 
alternating fields with pulse lengths of r ^ 0.5 — 0.6m~^. As a matter of fact, this is 
what one would expect from an investigation of the particle number of a single pulse 
as shown in Fig. 4. 22. The particle number shows a maximum at r ~ 0.5m~^ [28] 
corresponding to the threshold for creating the rest mass of the electron-positron 
system. Due to the fact that we have not chosen any constraint with respect to 
the oscillation frequency, it seems that the optimization procedure drives the field 
configuration towards the multiphoton optimum. 

In order to find a maximum of the particle number for any initial field config- 
uration, it seems to be best to add short pulses in alternating configuration with 
r ~ 0.5m^^. Hence, the search for an optimal field configuration can be reduced to 
the two cases presented in Fig.4.23. It can be seen that either the peak or the valley 
of the optimized short pulse configuration is at t = 0. Nevertheless, one has to be 
aware of the fact that the short-pulsed substructure will dominate the production 
rate and thus, mainly the multiphoton effect will contribute to the particle number. 
Accordingly, the optimal field configuration with regard to particle production does 
not correspond to the dynamically assisted Schwinger effect, but to the multiphoton 
effect. This is in agreement with previous studies, where multiphoton pair produc- 
tion becomes perturbative for r ~ 0.5m~^ [26] whereas the dynamically assisted 
Schwinger effect is still exponentially suppressed [32]. 




Time t[l/m] 

Figure 4.19: Deformation of an optimized electric field for time-dependent con- 
straints. The upper limit in the field strength is obtained at two points, whereas 
the lower limit does not affect the result. The initial field parameters are e = 
0.01, r = 5m~^ whereas the bounds are r = 4m~^ and r = 6m~^ and the offset 
Ae = ±0.0001. 
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Time t[l/m] 

Figure 4.20: Deformation of an initial electric field without constraints. The initial 
field parameters are e = 0.008, r = 5m~^. 




Figure 4.21: Deformation of an optimized field for different initiaV configurations 
and a constant constraint. One obtains similarities in the structure of optimized fields 
even for different initial fields. The initial configurations are e = 0.01, r = 3m~^ in 
a and e = 0.1, r = 20m~^ in c . In 6 and d the corresponding substructures are 
shown. Most notably, the oscillation rate and the peak height of the substructure 
resemble each other. 
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Figure 4.22: Dependence of the particle number on the pulse length for a single 
pulse in the multiphoton regime. For any given field strength the particle number 
increases with decreasing r until it reaches r ~ 0.5m~^, corresponding to the 
threshold for creating the rest mass of the an electron and a positron. 
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Figure 4.23: The two possible cases for a local maximum in the particle number 
within our approach. The oscillations are on top of the initial field with either 
a positive or a negative sign. The global maximum in the particle number then 
depends on the chosen constraints. 
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5 

Conclusions 



In this thesis we have presented an investigation of particle pair production in strong 
electric fields within the framework of quantum kinetic theory. We focused on field 
configurations consisting of superpositions of Sauter pulses. Our main findings are 
the following: 

We further investigated the dynamically assisted Schwinger effect. In compari- 
son to previous work, we focused on the effect of superimposing a strong and long 
pulse with a different number of short and weak pulses. In this respect, we have also 
examined the idea of a similarity between Ramsey-interferometry and the Schwinger 
effect. Our simulations revealed that the particle distribution function is sensitive 
to the specific form of the electric field. A further result of this analysis suggests 
that an indirect measurement of the Schwinger effect in well-adjusted experimental 
setups might be possible. 

In the course of our simulations it turned out that the combination of different 
solvers might enhance the numerical performance significantly. In this respect, we 
want to emphasize that a thorough selection of both solution and interpolation al- 
gorithms can save a large amount of computer time. 

Most notably we have derived equations to do pulse-shaping within the frame- 
work of optimal control theory systematically. The results of this investigation have 
revealed that the optimization procedure drives the system into a well-defined di- 
rection. Starting from a given initial field the system evolves into a state, in which 
the initial field is superimposed by an oscillating field. The oscillation rate of the 
superimposed field seems to correspond to the perturbative threshold for particle 
production. 
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This research has raised questions, which need to be further investigated. Unfor- 
tunately, an explanation for the deformations in the particle distribution as observed 
in our simulations cannot be given at the moment. Further studies on this subject, 
probably in complementary approaches, could reveal this puzzle. Moreover, further 
research might explore configurations with more realistic parameters in the optical 
or X-ray regime. Another focus should be laid on the inclusion of magnetic field 
effects and spatial inhomogeneities. Finally, we think that the idea of optimization 
to enhance pair production should be pursued in the future. Further investigations 
could offer new opportunities in this as well as in related research areas. 

All in all one can conclude there are still unsolved mysteries hidden in the area of 
strong field QED and, particularly, pair production. The rapid advances in LASER 
technology give us hope that an experimental observation based on explicit theoret- 
ical predictions comes into reach within the next decade. 
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Appendix A 

Additional Derivations 



A.l Normalization of Mode Functions 

The mode functions x"*" {q, t) and x~ {q^ t) form the fundamental system of the dif- 
ferential equation 

{d^ + co' {q, t) + ieE (t)) x {q, t) = 0. (A.l) 

For a vanishing electric field E (t) ^ the term u {q, t) becomes uj{q) and we obtain 
the simpler oscillator equations 

{d^+uj^q))xo{q,t) = 0. (A.2) 

The linear independent solutions of this equation are clearly given by plane waves 
and thus the mode functions are given by 

Xo^ iQ, t) = ^e-'-(^"^*, xo t) = ^e^'^(^^*, (A.3) 

where N"^ are normalization factors. They can be found by adopting the orthogo- 
nality relations for the spinors Ur{q,t) and Vr{—q,t) 

ul {q, t) ■ Us (g, t) = vl {-q, t) ■ Vs {-q, t) = Srs, (A.4) 
Ur {q, t) ■ Us {q, t) = -Vr {-q, t) ■ {-q, t) ^ —5rs, (A.5) 
ul{q,t)-Vs{-q,t)^0, (A.6) 

where r = 1,2 s = 1,2. For the sake of a simple presentation we just show the 
calculation of u\ ■ ui in explicit form here. The spinor Ur was defined by 

Ur{q,t)^{i^''dt-^-7r{q,t)+m)x^{q,t)Rr r = l,2 (A.7) 
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so we can write for vanishing electric fields 

Ur {q, t) = {h^dt - 7 • TT (g) + m) {q, t) Rr- 
Using the definitions for the gamma matrices we obtain 



(A.8) 



7 • TT 



17 



/ 


TTs TTi — m2 

\Ki + i7r2 -TTs 

/O i 0\ 

i 

1 
\0 i 0/ 



-TTa 
-TTi - 'm2 





-TTi + i7r2^ 








(A.9) 



(A.IO) 



Accordingly, the spinor takes the form 

Ui {q, t) 



( 



\ 



u 



m 


idt - q3 

\-qi -mj 

\ {q, t) = {m -idt - q3 -qi + 1^2) (x^)* {q, t) 



(A.ll) 



(A.12) 



Using the orthogonality relation between the spinors u^u — 1, their scalar product 
gives 



u 



[■ui = m'^ |Xo f - {idt + ^3) (Xo )* {idt - qs) Xo 

-{-qi + k2){xty {qi + iq2)xo (a.i3) 



and using the relations 

{qs + idt) {xoY {qs - idt) xt 

= ql IXol' - (dtxt) {Xty + iqsXt {dt (Xo+)*) + I^^Xol' (A-M) 

- {qi + m) xt {-qi + iq^) {xoY = ?i |xo f + ?2 |xo f 

one ends up with 



(A.15) 



u\-ui = |xo r + 9^ Ixc^r + 193 (xo {dt (xo)*)) 

-{{dtxt){xtT) + \dtxt\"- (A.16) 
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In combination with equation (A. 3) this leads to an equation for the normahzation 
constant N'^ 



-iijjt 



( icue'' 



7V+ ) (N+y 



+ 



-iojt 



N+ 



This is equal to 



and yields to the result 



u"^ — q^u — q^u + cu^ = \N~^ \ 



N+ (q) = y/2uj(q)(uj(q)-qs). 



1. (A.17) 
(A. 18) 
(A.19) 



Using the orthogonality relation between the spinors v'^v = 1, one obtains the nor- 
malization constant N~ — \/2u {oo + q^). Thus we can express both mode functions 
for a vanishing electric field as 



^2a; (q) {u (q) TQs) 



(A.20) 



It can be checked explicitly that all other orthogonality relations hold when we use 
this normalization. 

A. 2 Bogoliubov Coefficients 

We start with the spinor ip {q, t) in the particle and in the quasi-particle represen- 
tation 



r=l 
2 

= J2 {Ur (q, t) Ar {q, t) + Vr {-q, t) Bl {-q, t)) 



r=l 



By inserting the expressions (2.53a) and (2.53b) 



Ar {q, t) = a {q, t) a, (q) - 13* {q, t) bl (-g) , 
Bl i-q, t)=f3 (q, t) ar (q) + a* (q, t) bl (-q) 



(A.21) 
(A.22) 



(A.23a) 
(A.23b) 
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into the equation above one gets the expression 

2 

J2 ((7°^ - 7 • TT + m) K+Rr {aur - P*bl) 

r=l 

+ (-7°a; - 7 • 7f + m) K'Rr (/3a^ + a*bl)) 

2 

= X] ((i7°^t - 7 ■ + m) x^Rvar 
+ (i7°9t-7-7f + m) 



(A.24) 



We spht this equation now in different parts for better readabihty. The left hand 
side gives for r = 1 





\-7r1 - iTTa/ 



K ' I aai 



+ 





\-7r1 - iTTa/ 



(|/3ai + a*h^^ . (A.25) 



Consequently we find for r = 2 
/tti - i7r2\ 



a; - TTs 




The right hand side gives for r — 1 

f rUe \ 


idt - TTs 

\-7r1 -iTTa/ 



/ TTi - i7r2\ 




K- [pa2 + q;*4) ■ (A.26) 



/ rUe \ 


i^t - TTs 

\-7r1 - i7r2/ 



X 



(A.27) 



and for r = 2 



/tti - i7r2\ 


\ rUe J 



X 02 + 



/tti - iVTsX 
i^t - TTs 





X-bl 



(A.28) 



The obtained expression is a vector equation with four components. To keep the 
treatment simple, we do not deal with all four components at the same time. Rather 



76 



we start with the first component 

+ rrieK 

+ (tti - iTTs) (aa2 - /3*bl^ + (tti - iTTs) K- (^/3a2 + a*bl^ (A.29) 

= meX+ai + meX~b\ + (tti - ivra) x+Oa + (tti - ivTs) x"6^ 



X+, me 7^0, (A.30) 

X+, (7ri-i7r2)7^0, (A.31) 

X", me 7^0, (A.32) 

X", (7ri-i7r2)7^0. (A.33) 

We are deahng with massive particles so we can use equations (A.30) (A.32) and 
neglect the restrictions of the equations (A.31) (A.33). We proceed with the second 
component 

^{idt-7rs)x^a2 + {idt-7rs)x-bl (A.34) 

and again we find via comparison of coefficients 

02 : {u - TTs) K+a - {u + tts) P = {idt - tts) x"^ (A.35) 

bl: -{u;-tt3)k-^I3* -{cu + 7Ts)K-a* ={idt-ns)x-. (A.36) 

We get four different equations (A.30), (A.32), (A.35) and (A.36), as the lower two 
components in (A. 24) do not contain additional information. Now we restrict our- 
selves to the coefficients of the operators Oj and pick out the two equations 

K+a + K~I3 =x"^, (A.37a) 
{ou - TTa) K+a -{uj + TTs) K-p = {idt - TTg) X^- (A.37b) 

Multiplying equation (A. 37a) with (a; + tts) and taking the sum of both equations 
gives 

{uj + TTa + uj - TTs) K'^a = (a; + tts) x"^ + (i^* - tts) x"^, (A.38) 

which yields 

2ujK'^a = ux^ + i (9tx^) ■ (A.39) 



and find via comparison of coefficients 



oi : 


K'^a + K 




a2 : 


K'^a + K 






-K+P* + K- 


a* 




-K+P* + K- 


a* 
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Then we multiply this equation by {2ujk'^) ^. As was defined as 

-i9(g,t) 

K+ {q, t) = . , (A.40a) 

^y2uJ {q, t) {U (g, t) - TTg {q, t)) 

ie{q,t) 

K- {q, t) = , , (A.40b) 



we can write 



,2 ^2 



Thus, we obtain for a {q, t) 

a {q, t) = i^^2 {q,t)-7:l{q,t)K- {q, t) {dt - icu {q, t)) x+ {q, t) . (A.43) 

A similar calculation holds for (3 {q,t), too. Again we start with the equations 
(A. 37a) and (A. 37b), but this time we multiply equation (A. 37a) with — (a; — tts). 
The resulting equation is then 

2UJK-/3 = {u -idt)x^- (A.44) 
After a calculation similar to the one above we obtain 



/3 (g, t) = -i^u^ {q, t) - 7r| (g, t)K+ (g, t) {dt + icu (g, t)) x+ (g, t) . (A.45) 
A. 3 Derivatives of the Bogoliubov Coefficients 

Wc start with the definition of a {q, t) as given in (A.43) but explicitly insert the 
expression for k~ 



a 



{q, t) = i ((a,x+ {q, t)) - ico {q, t) x+ {q, t)) (A-46) 
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Then we differentiate a {q, t) with respect to t 



+ i((9a-^)-ia;x+)a,|^/^e^^). (A.47) 



(*) 

The time derivative of the term (*) gives 



1 / 2uj (d; - Tfa) 2a; - (a; - TTs) 2a; -g uj - rcs -.g. 

e + \ — e loj. 



(A.48) 

2\l oj -TTs \.bP- V 2a; 

Taking into account the expressions 

a; = (A.49) 
Tfs = eE, (A.50) 

we can further rewrite the term (*) as 



1 / 2a; eE-Ks — eEu — eE7rs + eEnl/u u — tts -^g. 

-e + \ — e la; 



2\uj — 7rs 2uj^ V 2a; 



(A.51) 

Inserting (A.51) into (A.47) we obtain upon taking into account the differential 
equation (A.l) 



= ly e ( -a; X - le^^X - ^^^^^ ~ ) 

-(.x-i.x^)(i.^.i.)). 



(A.52) 
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This can be simplified to 



• 1^ le [ [ • c -e-ETTs \eE ieETTsN + , 



1 



(A.53) 



with (g, = ^51^ and el (q) ^ml + ql + ql 

Analogously, we perform the same lengthy but straightforward calculation for 



' t) = -i (g, t)) + ia; t) x+ t)) W ^(g>^)+J^3(g,t) ^_,,(,,)_ 



2a; 



Differentiation with respect to t gives 



- i + ic.X+) a, (/^e-^] . (A.55) 



(**) 

The time derivative of the term (**) gives 



^ • _ (A.56) 

^ -10 / _■ \ 

2 V + TTa 4c<j2 ^ ' V 2a; ~^ ^ ^'^'^ 

Taking into account the expressions uj — ^^^^^ — eE again, we can further rewrite 
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the term (**) as 

1 I 2uj cEtts + eEui — eETTs — eEir'^/uj Iu + tt^ 



2 V + TTs 2cc;2 V 2w 



We put everything together and get 



(jj + Tiz _i0 f f . ^ ^ . ^uj , . e£' . _ eE\ 



+ ((^a-)g(--3))) 



= --W^ae-2i^. 
2 

This completes the check on the Bogohubov coefficients. 



(A.57) 



2a; VV TTs 2 20;/ (^.58) 



A. 4 Hypergeometric Differential Equation 

We start with the differential equation 

(^u (1 - u) duU (1 -u)dy, + ^uj^ (g, u) + ier^u (1 - u)^ % (g, m) = (A.59) 
and the ansatz 

X{q,u)^u''{l-ufrj{q,u). (A.60) 

In order to solve the differential equation (A.59), we rewrite it by using the ansatz 
(A.60) starting with the derivatives. The first derivative with respect to u is given 

by 

(duX) = (aw"-' (1 - uf - f3u'' (1 - uf-^^ X 77 

+ (u''il-uf^ x{duv) (A.61) 
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and the second one is given by 

+/3u'' {(5 - 1) (1 - uf-^^ X ri 
+2 [au'^-^ (1 - uf - pu'^ (1 - uf-^^ X [d^T]) 

+ (1 - uf^ X a„ (duV) . (A.62) 

Note that we have to apply the product rule in the first term in (A. 59) 

u{l-u) duU (1 - u) (duX) 

^u{l- uf (a„x) - (1 - u) (duX) + uf (a„ {duX)) ■ (A.63) 

We then insert these derivatives into the differential equation (A. 59) resulting in a 
lengthy expression. The prefactor of the second derivative du {dyrj) is given by 

ii"+2(l-M)^+', (A.64) 

whereas the prefactor of the first derivative [dyj] {q,u)) reads 

2a«"+^ (1 - uf^' - 2/3«"+2 (1 _ 

+ (1 - uf+^ - (1 - uf+^ . (A.65) 

Finally, the prefactor oir]{q, u) is the most comphcated one. There are three contri- 
butions in total, one arising from (A.62) 

a (a - 1) (1 - uf+^ + /3{/3-l) (1 - uf - 2a/3M"+^ (1 - uf^' , (A.66) 
another one stemming from (A.63) 

au"" (1 - uf-^^ - (1 - uf+^ - (1 - uf+^ + (i _ ^f (A.67) 

and a last one originating from the remaining terms in (A. 59) 

2 

^ (a;2 + 4isu (1 - u)) u" (1 - uf 

^^{l + ql + ql + ql + 2^3^™ - 2q^eT (!-«)+ ffV^ ^^"^^^ 
+£V (1 - uf - 2£V'U (!-«) + to (1 - uf (1 - uf , 
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where we have used 

{2u - if = Au^ - Au + I = - 2u{l - u) + {I - uf . (A.69) 
Accordingly, the differential equation eventually is given by 

«"+2 (1 _ „)/3+2 X 3^ (a„r; (g, u)) (A.70) 
+ (^2aM"+i (1 - uf^^ - 2/3u''+^ (1 - uf^^ 

+u^+' (1 - uf"-^ - u'^^^ (1 - uf-^^') X (9„77 (q, u)) 
+ (a (a-l)u''(l-uf^^ 



+/3 {/3 - 1) u'^+^ (1 - - 2a/5M"+i (1 - uf+^ 
+au'' (1 - uf+^ - /3u''+^ (1 - uf^^ 

+— {l + ql + ql + ql + 2q^eTU - 2q^eT {I - u) + s^t'^u 

+£^2 (1 - uf - 2e^T\ (1 - ii) + 4i£'u (1 - u)) u'^ (1 - uf x r] {q, u) (A.71) 

= 0. 

A. 5 Normalization of Analytic Mode Functions 

The normalization constants are chosen such that 

(u) ^K^(u), u^ 0+. (A.72) 

Accordingly we have to investigate the behaviour of the mode functions x {Qj in 
this limit 

X+ {q, u) = N+u"" {1-uf 2^-1 (a, b, c; u) 

as 2^1 {a,b,c;u 0+) = 1. A similar calculation can be done for x~ {q^u), which 
results in 

X~ (g, u) = N-u-'^, u 0+. (A.74) 
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It is possible to calculate k {q, u) in the same limit. Prom the definition of the 
functions, we have 

g=Fi6»(g,w) 

(g, u) = , (A. 75) 

^y2u (g, u) {u {q, u) =F 713 (g, u)) 

with the transformed dynamical phase 9 {q, u). Using dt — ^u{l — u) we find 

9{q,t)^ [ uj{q,r)dT, (A.76) 

J to 

which takes the form 

^(,»=2 r^^'i^^. (A.77) 



^Juo U{l-U)' 

The integration can be performed analytically 

— — — ^ rdu — — 2£Tln (^3 + STU — ST {1 — u) (g, u)) 

+ Wo (In {u) - \Yi{l + ql + ql + {q^ - er) (gg + (2m - 1) er) + UqU (g, u))) 

- ui (In (1 - m) - In (1 + g? + ga^ + (ga + £r) (gs + (2m - 1) st) + c^io; (g, u))) , 

(A.78) 

where uq = oj (g, 0) , cji = w (g, 1). In the limit u — ;> 0+, the dynamical phase splits 
into a regular part 9 (g, 0) as well as into a divergent part 

9 {q, u) ":i>°^ ^a;oln (u) + 9 {q, 0) . (A.79) 

Denoting 773 (g, 0) = g3 — st the mode functions k read 

(q, 0) = ^^T^^oHu)_ (A.80) 

^200 {q,0) {oo{q,0)T7rs{q,0)) 

Comparison of (A. 73), (A. 74) and (A. 80) fixes the normahzation constants 

N+ (g) = , (A.81) 
^2u{q,0) {u (g',0)-7r3(g',0)) 

N~ (q) = , (A.82) 

^2u{q,0) {u (g,0) + 7r3 (g^O)) 



84 



A. 6 Asymptotic Distribution Function 

We found an analytical expression for the distribution function 

F%u) = |iV+(g-)f fl + ^^iigl^) \F,iq,u)+iF,iq,u)\\ (A.83) 

(*) 

where 

Fi (g, u) = -u (1 - u) — 2 -Fi (1 + a, 1 + 6, 1 + c; m) , (A.84) 
r c 

F2 {q, u) = {00 {q, u) - {1- u) ooq - uuji) 2^1 (a, c; u) . (A.85) 

We now want to derive an asymptotic expression for the distribution function in the 
limit M — > 1~. As we will encounter the gamma function r(a;) in the following, we 
introduce it here. This gamma function is defined for all complex numbers except 
the non-positive integers. For n G N the gamma function coincides with the factorial 
function 

r(n) = (n-l)!. (A.86) 
The gamma function obeys, among others, the following relations 

r(i + .) = .rw, Ir(i + i.)r = ^. (a.st) 

In order to determine the limit m — )■ 1^ of the hypergeometric function 2^1 (a, 6, c; ^), 
we may transform it according to linear transformation rules [54] 

^ , , , r(c)r(c-a-6) ^ , , 

27^1 (a, 6, c; z) = — — — —27^1 (a, 6, a + 6 - c + 1; 1 - 

V [c — a)V [c — b) 

+ (1 - zf""'^ ^^^^V{a)V^) ^^ ^"^^ (c - a, c - 6, c - a - 6 + 1; 1 - ^) . 

|arg(l-z)| < TT (A.88) 

We will focus on the term (*) in (A.83) for the moment as it is the most compli- 
cated one. In the limit — > 1~ we obtain 

Fi (g, u) = (1 - u) — 2^^! (1 + a, 1 6, 1 c; m) 
T c 

2 a& c-a-6-1 r (1 + c) r (1 + a -I- 6 - c) 
= -M 1 — M — 1 — M — ^ — \ — -, TT A. 89 

T ^ ^ c ^ ' r(i + a)r(i + 6) ^ ' 
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and 



F2 {q, u) = {uj {q, u) - {1 - u) uq — vwi) 2J^i {a, b, c; u) 

w^i- //I N /I \ \l^ vc-a-6 r (c) r (a + 6 - c) 
= {{1-u)ujx-{1-u)ujq){1-u) p ^ ■ (A.90) 



c—a—b 



This term 



In the equation |Fi + iF2|^ we can single out the term {1 — u) 

gives always 1, even in the limit w — )■ 1~. The second term in the equation above 
vanishes completely in this limit (because of the remaining {1 — u) term) leading to 



|Fi + iF2 



2 u->l" 



abr(l + c)r(l + a + b-c) 



c r(i + a)r(i + 6) 



(A.91) 



This then yields 



«->-i 



oJi + TTs {q, 1) 4 



2uo{uo - TTsiq^O)) Ui 



abr{l + c)T{l + a + b-c) 



c r(i + a)r(i + 6) 



(A.92) 

We define the following shorthand notations for the arguments of the gamma func- 
tions 



1 + a = 1 — 1 er H — — - 



+ + + (A.93) 



(A.94) 



C— 1 — ITUJo = 1+1 {—TUJo) = 1 + 17 



(A.95) 



1 + a + b- c^l + irooi = 1 + i5. (A.96) 
In a next step we rewrite the asymptotic distribution function according to 



F{q,u^l ) = 



(^1 +71"! 2 

a 



r (c) r (1 + a + 6 - c) 



r (1 + a) r (6) 



2 + TTi 2 vr7(5sinh (ttq;) sinh (7r;5) 

(X 



t'^cUqCuI cuq — TTo TT^CK^sinh (777) sinh (7r(5) 



(A.97) 
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F{q,u^l ) = p 

— 2 p 

LOq — TTo 2eT — Uq + Ui 



where we have used 



UqCOi — TTqUJi + TTiUJo — 2eTT^Q — Uq 



TTo = g3 - ST, 771=^3 + £T, (A.99) 

ul + TT^, + TT^, (A.lOO) 

(A.lOl) 



and 



_ sinh (ttt {2eT + - wi) /2) sinh (ttt {2eT - o-o + a-i) /2) 
^ sinh (tttcuq) sinh (7rTa;i) 

We further simphfy the following expressions 

2£r7ri -Luj = -£]_ - ql- 2eTqz - (er)'^ + 2eT {qs + er) = 

-el-ql + ierf, (A.103) 



- 2£r7ro - = -si - ^3 + 2£t?3 - {erf - 2eT {qs - er) = 

-el-ql + ierf, (A.104) 

such that 

25r7ri - cjj = -2£r7ro - a;o. (A. 105) 

Altogether, we obtain the final result 

j^u^i- r, _ sinh (ttt (2£T + Wo - /2) sinh (vrr (2£r - a;o + /2) /..^^x 

^ Zp^z — — ^— r-r-7 ^ ■ (A.iUbj 

smh (7rTa;oj smh [tttuji) 
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Appendix B 

Runtime Tables 



To provide additional information concerning the numerical performance, we give 
here the unnormalized data from section 3 on Numerics. We have done all test 
calculations on an "Intel(R) Core(TM)2 Duo CPU E7400 @ 2.80GHz" with the 
GNU-C++-Compiler g++-4.5. For testing we tried to use always two threads at the 
same time. The templates for the C++ code were obtained from the books [55, 56] 

Table B.l: Runtime comparison for different field configurations and 
solvers (c.f. Tab. 3.1). 



Field Parameters 
ie[l], r[l/m], tO[l/m]) 


Standard 
(sec) 


Standard(I) 
(sec) 


Trans2(I) 
(sec) 


Trans2 
(sec) 


Transl 
(sec) 


£ = 0.1, r = 100^ 


292.54 


197.72 


201.62 


266.4 


647.11 


£ = 0.005, T = 3^ 


8.79 


6.56 


6 


2.75 


13.29 


£1 = 0.1, Tl = 100, 
£2 = 0.01, T2 = 22 


72.03 


65.63 


68.38 


112.67 


326.75 


£1 = 0.1, Tl = 100, 

£2 = 0.01, T2 =2,to = 30^ 


421.34 


313.09 


257.44 


376.56 


868.15 


£ = ±0.005, T = 3,to = 53 


343.44 


231.67 


216.52 




734.49 


£ = 0.005, T = 3, to = 53 


128.67 


108.08 


219.61 


187.66 


231.99 


£1 = 0.1, Tl = 100, 

£i = ±0.005, Ti = 3, to = 9* 


101.21 


83.73 


186.8 


179.1 


223.02 


£1 = 0.1, Tl = 100, 

£i = 0.005, Ti = 3,to = 9* 


830.95 


543.5 


350.96 


463.78 


1839.54 


£1 = 0.1, Tl = 100, To = 200, 
e, = 0.005, Ti = 3, to = 24.95"' 


829.96 


540.29 


357.01 


450.81 





^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long 
pulse and ten short pulses 
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Table B.2: Comparison of the particle number for different solvers and field con- 
figurations (c.f. Tab. 3. 2). 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


Standard 


Standard(I) 


Tr2(I) 


Trans2 


Transl 


An 






£ = 0.1, T = 100^ 


1.529 


1.668 


1.529 


1.540 


1.529 


1.529 


xE- 


13 


£ = 0.005, T = 3I 


5.856 


5.856 


5.856 


5.933 


5.856 


5.856 


xE- 


10 


£ = 0.05, r = 100^ 


1.392 


9856. 


1.193 


0.599 


1.392 


2.4E-09S 


xE- 


18 


£1 = 0.1, Tl = 100, 

£2 = 0.01, T2 = 2^ 


4.240 


4.240 


4.240 


4.240 


4.240 




xE - 


07 


£1 = 0.1, Tl = 100, 

£2 = 0.01, r2 = 2, to = 30^ 


4.021 


4.021 


4.021 




4.021 




xE — 


07 


e = ±0.005, T = 3, to = 5'^ 


8.700 


8.700 


30.99 


31.03 


8.700 




xE - 


10 


e = 0.005, T = 3, to = 53 


3.782 


3.782 


5.241 


5.267 


3.782 




xE- 


10 


£1 = 0.1, Tl = 100, 

£i = ±0.005, Ti = 3,to = 9"* 


1.421 


1.421 


1.421 


1.420 


1.421 




xE - 


08 


£1 = 0.1, Tl = 100, 

£i = 0.005, Ti = 3, to = 9" 


1.353 


1.353 


1.353 


1.353 






xE- 


08 


£1 = 0.1, Tl = 100, To = 200, 
£1 = 0.005, n = 3, to = 24.95'' 


2.259 


2.259 


2.259 








xE- 


08 



^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long pulse and ten 
short pulses ^ Discrepancy coming from underflow errors 



Table B.3: Runtime comparison for different mo- 
mentum spacings for a single pulse with e — 
0.005, T = 3m~^. All momenta are given in units 
of the electron mass q[m] (c.f. Tab. 3. 3). 





Aqx = 0.5 

(sec) 


Aq± e [0.1,0.5] 
(sec) 


(sec) 


Important qs ^ 


0.95 


1.55 


2.91 


Aqs = 0.05 


1.15 


1.54 


3.10 


Aqs = 0.01 


3.62 


6.62 


59.57 



^ Ag3 = 0.05, 0.1 Tolerance 
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Table B.4: Runtime comparison for different mo- 
mentum spacings for a single pulse with e — 
0.1, T = lOOrn"^. All momenta are given in units 



of the electron 


mass g[m] 


(c.f. Tab. 3.4). 






Agx = 0.5 
(sec) 


Agx e [0.1,0.5] 
(sec) 


Agx = Aga 
(sec) 


Important qz ^ 


189.43 


403.47 


655.78 


Ags = 0.05 


254.82 


557.61 


917.52 


Aij3 = 0.01 


1181.22 


2663.16 


18685.5 


^ Agg = 0.05, 0.1 


Tolerance 







Table B.5: Runtime comparison for different sets of sampling points. In 
order to calculate intermediate values linear interpolation is used (c.f. 
Tab.3.6). 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


Standard 
(sec) 


Trans2(I) 
(sec) 


512 Points 
(sec) 


4096 Points 
(sec) 


32768 Points 
(sec) 


£ = 0.1, r = 100^ 


292.54 


201.62 


74.21 


229.4 


740.98 


£ = 0.005, r = 3I 


8.79 


6 


2.91 


39.3 


201.14 


£1 = 0.1, Tl = 100, 
£2=0.01,r2=22 


421.34 


257.44 


85.41 


222.17 


800.55 


£1 = 0.1, Tl = 100, 

£2 = 0.01, T2 = 2,io = 30^ 


343.44 


216.52 


75.65 


180.82 


643.2 


£ = ±0.005, T = 3, to = 53 


128.67 




12.8 


180.94 


571.54 


£ = 0.005, r = 3, to = 53 


101.21 




5.19 


155.36 


507.67 


£1 = 0.1, n = 100, 

£i = ±0.005, Ti = 3,to = 94 


830.95 


350.96 


75.2 


138.99 


581.54 


£1 = 0.1, Tl = 100. 

e, = 0.00.')^ T, = :-!. /„ = !)-^ 


829.96 


357.01 


71.93 


119 


575 



^ Single; pnlsc' ^ Double pulse ^ Ten identical pulses with time lag ^ One long pulse 
and ten short pulses 
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Table B.6: Comparison of the particle number for different sets of sampling points 
In order to calculate intermediate values linear interpolation is used (c.f. Tab. 3. 7). 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


Standard 


Trans2(I) 


512 Points 


4096 Points 


32768 Points 






£ = 0.1, r = lOQi 


1.529 


1.529 




1.58659 


1.52987 


xE 


- 13 


£ = 0.005, T = 3I 


5.856 


5.856 


5.73369 


5.85943 


5.85716 


xE 


- 10 


£1 = 0.1, Tl = 100, 
£2 = 0.01, T2 = 22 


4.240 


4.240 


2.21893 


2.74849 


4.21281 


xE 


— 07 


£1 = 0.1, Tl = 100, 

£2 = 0.01, T2 = 2,to = 302 


4.021 


4.021 


1.67201 


2.57259 


3.99372 


xE 


-07 


£ = ±0.005, T = 3, to = 53 


8.700 




8.01398 


8.69707 


8.70156 


xE 


- 10 


£ = 0.005, T = 3, to = 53 


3.782 




3.43572 


3.78036 


3.78314 


xE 


- 10 


£1 = 0.1, ri = 100, 

£i = ±0.005, Ti = 3, to = 94 


1.421 


1.421 




1.00616 


1.41359 


xE 


-08 


£1 = 0.1, Tl = 100, 

£i = 0.005, Ti = 3,to = 9* 


1.353 


1.353 


7.92583 


9.19399 


1.34551 


xE 


-08 



^ Single pulse ^ Double pulse * Ten identical pulses with time lag to One long pulse and ten 
short pulses 



Table B.7: Runtime comparison for different sets of sampling points. In 
order to calculate intermediate values cubic interpolation has been used 
(c.f. Tab.3.8). 



Field Parameters 
(£[1], T[l/m], t0[l/m]) 




Standard 
(sec) 


Trans2(I) 
(sec) 


512 Points 
(sec) 


4096 Points 
(sec) 


32768 Points 
(sec) 


£ = 0.1, T = 100^ 




292.54 


201.62 


282.27 


283.29 


325.25 


£ = 0.005, T = 3I 




8.79 


6 


10.4 


13.94 


50.41 


£1 = 0.1, Tl = 100, 
£2 = 0.01, T2 = 22 




421.34 


257.44 


296.55 


307.23 


401.08 


£1 = 0.1, Tl = 100. 
£2 = 0.(11. r2 = 2.^0 = 


:-!()- 


34.'^. 44 


216.52 


236.8 


249.4 


312.73 


£ = ±0.005, T = 3, to = 


53 


128.67 




25.85 


29.84 


58.84 


£ = 0.005, T = 3,to = 5 


3 


101.21 




18.47 


24.25 


56.42 


£1 = 0.1, Tl = 100, 

£i = ±0.005, n = 3, to 


= 94 


830.95 


350.96 


188.4 


256.49 


271.99 


£1 = 0.1, Tl = 100, 

£i = 0.005, Ti = 3, to = 


94 


829.96 


357.01 


170.44 


255.34 


274.85 



^ Single pulse ^ Double pulse ^ Ten identical pulses with time lag to ^ One long pulse 
and ten short pulses 
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Table B.8: Comparison of the particle number density for different sets of sampling 
points. In order to calculate intermediate values cubic interpolation has been used 
(c.f. Tab.3.9). 



Field Parameters 
(£[1], r[l/m], tO[l/m]) 


Standard 


Trans2(I) 


512 Points 


4096 Points 


32768 Points 






£ = 0.1, r = 100^ 


1.529 


1.529 


1.4597 


1.51973 


1.52846 


xB - 


13 


£ — U.UUO5T — 


O.OOD 


O.oOD 


0.01 / Dl 


u.c5oCS41:4 


O.o04o 


Xr!/ — 


1 n 


£1 = 0.1, ri = 100, 

£2 = 0.01, T2 = 22 


4.240 


4.240 




3.16232 


4.22142 


xE- 


07 


£1 = 0.1, Tl = 100, 

£2 = 0.01, T2 = 2,to = 30^ 


4.021 


4.021 




2.97061 


4.0021 


xE - 


07 


£ = ±0.005, r = 3,io = 53 


8.700 




7.99011 


8.63437 


8.69277 


xE - 


10 


£ = 0.005, r = 3, to = 53 


3.782 




3.63728 


3.79134 


3.78425 


xE - 


10 


£1 = 0.1, n = 100, 

£i = ±0.005, Ti = 3, to = 9* 


1.421 


1.421 




1.12703 


1.41548 


xE - 


08 


£1 = 0.1, n = 100, 

£i = 0.005, Ti = 3,to = 9'' 


1.353 


1.353 




1.0406 


1.34787 


xE - 


08 



^ Single pulse ^ Double pulse * Ten identical pulses with time lag ^ One long pulse and ten 
short pulses 
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